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This Annual Status Report, Volume I, describes the results of work performed 
on Task VB (Special Finite Element Models) during the fourth year of the NASA 
Hot Section Technology (HOST) program, "3-D Inelastic Analysis Methods for Hot 
Section Components" (Contract NAS3-23697). The finite element code, MHOST 
(MARC HOST), is developed further based on mixed iterative solution procedures 
wFose concepts and basic technology were established in the first and second 
year efforts. The technology is extended in the current phase to incorporate 
the effect of multiple embedded singularities in generic modeling regions. 
Specifically, the local mesh refinement technology and the special function 
representation for singular behavior of strain are developed based on mixed 
element concepts. The local mesh refinement algorithm, referred to as the 
subelement iteration, utilizes conventional linear and higher order polynomial 
representations for spatial discretization. 

The program is being conducted under the direction of Dr. C. C. Chamis of the 
NASA-Lewis Research Center. Prime contractor activities at United Technologies 
Corporation are managed by Dr. E. S. Todd. Subcontractor efforts on finite 
element tasks at MARC Analysis Research are led by Dr. J. C. Nagtegaal . 
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1 . 0 SUMMARY 


The Special Finite Element Models portion of the HOST (Hot Section Technology) 
3-D Inelastic Analysis Methods program is divided into two 24-month segments: 
a base program, and an option program exercised at the discretion of the 
Government. Versions 1 and 2 of MHOST (MARC- HOST ) were developed during the 
base program (Tasks IB and I IB). The MHOST code employs both shell and solid 
(hexahedral) elements in a mixed iterative solution framework. This code 
provides comprehensive capabilities for investigating local (stress/strain) 
and global (vibration and buckling modes) behavior of turbine engine hot 
section components. In the development of MHOST, advantage has been taken of 
the technical expertise of engineers at MARC Analysis Research Corporation. 
This has led to the construction of sophisticated algorithms and codes which 
may reduce computer time and memory requirements for three-dimensional 
analyses. 

Version 4 of the MHOST code has been developed under Task VB of the HOST 3-D 
Inelastic Analysis Methods option program. These latest development efforts 
have concentrated on the consolidation of the numerical technology created and 
implemented for the basic mixed iterative solution strategy. In addition, the 
subelement solution strategy designed to capture local behavior related to 
embedded discontinuities has been investigated further. 

Efforts have been devoted to maintain the computer program on different 
computer systems and to improve quality. An in-core profile solution subsystem 
has been added, replacing the original band matrix solution. The computational 
procedures for eigenvalue extraction and transient time integration have been 
improved and been made compatible with this new solution package. 

Enhancement of the solution capabilities of the MHOST code include the large 
displacement option to update nodal coordinates, the stress stiffening option 
for quasi -static and eigenvalue analysis and the centrifugal mass option for a 
structure rotating at a high angular velocity. 

Algorithms developed for and implemented into the MHOST code utilizing the 
mixed global solution strategy and subelement techniques are discussed in 
detail in this report in conjunction with remarks on programming 
considerations. 


2.0 LITERATURE SURVEY 


This section reviews literature on the application of nonlinear finite element 
methods to the inelastic analysis of turbine engine hot section components. It 
is not intended to discuss here the state-of-the-art finite element technology 
but to provide a snapshot of the current trends in research and development 
for the period since the last annual report [NAKAZAWA (1986)]. An attempt is 
made in this survey to position the technology developed under the HOST 
contract in the computational mechanics perspective. 

2.1 INTRODUCTION 

The major areas of research in the improvement of computational performance 
and applicability of finite elements have continued to be: 

1. Mixed and hybrid formulations designed to improve the accuracy of 
finite elements 

2. Global solution procedures for linear and nonlinear finite element 
equations to reduce both memory requirements and the number of 
arithmetic operations with particular emphasis on iterative solution 
strategies 

3. Post-processing methods designed not only to estimate errors 
accurately a posteriori, leading to adaptive mesh refinement 
strategies, but also to extract more information out of a finite 
element solution. 

Both modern numerical technology and advanced computing machinery now allow 
the application of finite elements to a wider range of problem areas. The use 
of large meshes for representing fine detail in engineering problems has 
become very valuable. Sophisticated numerical methods have been developed to 
deal with combined nonlinear kinematics and material responses in a highly 
efficient manner. 

This section concentrates on the reports and articles which have appeared 
since the previous literature survey reports [NAKAZAWA (1986) Section 3.1, and 
FYHRE, HUGHES (1983)]. 

2.2 VARIATIONAL FORMULATION AND ELEMENT TECHNOLOGY 

Interest continues to remain high in using the mixed variational formulation 
as a basis for developing finite element methods. The Hu-Washizu principle has 
become the most popular variational formulation from which numerical methods 
are derived. It seems the most natural approach, for instance, from which to 
derive algorithms for probabilistic structural analysis. Evaluation methods 
for finite elements based upon mixed variational principles have also been 
developed recently. In modern literature, the words 'mixed' and 'hybrid', as 
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applied to finite elements, are often used interchangeably and ambiguously. 
Seemingly, the mixed methods derived from piecewise discontinuous stress 
representations are usually referred to as 'hybrid' formulations, while those 
derived from independent strain interpolations as well as other mixed 
formulations are referred to as 'mixed' forms. It is interesting to note 
moreover that the word 'hybrid' is even used by a certain sector of 
computational mechanics as referring to a class of methods combining the 
finite element and boundary integral formulations. 

Interest in the Hu-Washizu principle has been revived in recent years by the 
untimely death of one of the original investigators and the subsequent 
publication of his memorial volume [PIAN (1984)]. Utilization of this 
principle for the analysis of existing methodology, such as the numerically 
integrated displacement formulation [ZIENKIEWICZ, NAKAZAWA (1984)], appears 
frequently in the literature. BELYTSCHKO (1986) presents an excellent review 
of the technology available for thick plate and shell element formulations 
based on the Reissner-Mindlin theory. 

Efforts to use the mixed Hu-Washizu finite element equations as a driving 
mechanism for the iterative solution of linear and nonlinear problems continue 
as a part of the HOST project, as reported in progress reports by WILSON, BAK, 
NAKAZAWA, BANERJEE (1984, 1985), and NAKAZAWA (1986). Also, a number of papers 
discuss various aspects of this methodology [NAKAZAWA, NAGTEGAAL (1986), and 
NAKAZAWA, SPIEGEL (1986)]. 

With particular application to transient dynamics, a version of the mixed 
iterative solution algorithm combined with a second order, single step time 
integration operator has been proposed and applied to a number of linear 
elastic problems by ZIENKIEWICZ, LI, NAKAZAWA (1986). 

An analysis of a certain class of mixed finite elements is discussed by LE 
TALLAC (1986). An iterative algorithm based on the Hu-Washizu principle is 
derived for rate-independent deviatoric plasticity. The elastic strain is used 
as a variable instead of the total strain. This formal alteration of the 
formulation, in conjunction with the radial return algorithm being hardwired 
to the analysis, enabled Le Tallac to discuss the solution existence question; 
this then allows him to consider the possible formalization of convergence 
characteristics. 

For the evaluation of element formulations, as well as being the first step of 
finite element code validation, the patch test invented by Irons and 
documented in BAZELEY, CHEUNG, IRONS, ZIENKIEWICZ (1966) has been used by many 
researchers and commercial code developers. For instance, NAGTEGAAL, NAKAZAWA, 
TATEISHI (1986) discuss the development of an eight-node shell element based 
on the assumed strain approach and demonstrate its validity by performing an 
extensive set of patch test calculations. The concept of patch test is 
re-examined in a recent paper by TAYLOR, SIMO, ZIENKIEWICZ, CHAN (1986), where 
the equivalence between the satisfaction of the patch test and the consistency 
and stability of finite element approximations is discussed. Note that the 
argument by STUMMEL (1980) seriously questioning the validity of the patch 
test is countered in this paper. 
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An important extension of the patch test is proposed in a more recent paper by 
ZIENKIEWICZ, QU, TAYLOR, NAKAZAWA (1986). Assuming that the consistency of the 
approximations is satisfied, the test focuses upon the stability of mixed 
finite elements implied by the Babuska-Brezzi condition [BABUSKA (1973), and 
BREZZI (1974)]. The procedure involves counting the minimum possible number of 
active degrees -of -freedom in a patch and comparing this number with the 
maximum possible number of constraint degrees-of-freedom. If no constraints 
are imposed on nodal stresses and strains, the mixed iterative solution 
procedure with continuous interpolation passes the test and is assumed stable. 
An extension of this stability patch test to the three field formulations of 
Hu-Washizu type and to plates and shells derived from the Reissner-Mindlin 
theory is being investigated by ZIENKIEWICZ (1986) and his co-workers. 

The relation between the stability patch test and the mathematically derived 
Babuska-Brezzi condition is also being actively investigated. 

Use of the most general mixed variational formulation of the Hu-Washizu type 
in element development continues to be in the mainstream of research 
activities. Despite the slight confusion in terminology, a family of simple 
high performance elements has been developed utilizing, consciously or 
unconsciously, the mixed method concepts. 

Interest in the equal order interpolation for mixed finite element processes 
is definitely growing. A number of papers and reports have been published not 
only via this project but also by the research group at Swansea. Also, a 
number of associated academic research projects have been initiated in the 
last couple of years. It is anticipated that the number of publications on the 
utilization of the equal order interpolation mixed method and its iterative 
solution strategy will increase considerably in the next few years. 

2.3 A POSTERIORI ERROR ESTIMATES AND POST-PROCESSING 

The methodology used for computing the approximation error in a finite element 
solution has been improved significantly in the last few years. Also, the 
adaptive mesh refinement algorithms based on the spatial distribution of error 
indicators have been developed and applied to a wide range of linear and 
nonlinear problems. A book on this subject, compiled by Gago, represents 
present state-of-the-art technology [BABUSKA, ZIENKIEWICZ, GAGO, OLIVEIRA 
(1986)]. 

The papers included in the book indicate that a posterior error estimates and 
adaptive refinement methodologies are well understood mathematically and 
tested for a wide range of linear elastic problems. However, the computational 
processes look overly complicated and extensions to nonlinear analysis of 
solids and structures will need further research and development. Application 
of the adaptive mesh refinement concept has been demonstrated to be tractable 
using a relatively simple error indicator for solutions to aerospace flow 
problems modeled by the compressible Euler equations. Note that the adaptive 
solution of the compressible Euler equations is the first systematic attempt 
at extending the technology to nonlinear problems. The results included in the 
book show the potential advantage of the concept. 
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A simple algorithm to estimate errors a posteriori is proposed by ZIENKIEWICZ, 
ZHU (1987). The logic is based on the difference in values for quantities 
obtained directly within the elements and the same information projected by 
nodal interpolation functions. The nodal projection techniques used in the 
mixed, iterative process to recover the nodal strain allow a simple and 
effective way to calculate errors. This process is relatively easy to 
implement, in comparison with other seemingly complicated mathematically 
derived a posteriori error estimate algorithms. Examples shown in the paper 
indicate the efficiency of the proposed process. However, rigorous 
mathematical analysis is still to be done to verify the methodology. It is 
interesting to note that the residual vector calculated at the zeroth 
iteration of the equal order mixed iterative solution is indeed the 
Zienkiewicz-Zhu error indicator weight-averaged at nodes. 

With triangular elements being employed, efforts to combine automatic 
mesh-generation and adaptive refinement concepts have been reported. Examples 
included in ZIENKIEWICZ, ZHU (1987) use this idea. In the application of 
adaptive mesh refinement to compressible flow problems, a number of papers 
have been published in which triangular mesh generation methods are discussed 
in conjunction with local refinement strategy to capture shocks. 

PERAIRE, MORGAN, ZIENKIEWICZ (1986) and LOHNER (1986) describe adaptive mesh 
refinement and de-r^finement algorithms designed for triangles and speculate 
on the applicability of such a concept to three-dimensional computations. The 
use of automatic mesh generation of triangular and tetrahedral elements has 
been found to be viable in finite difference flow computations involving 
unstructured grids. JAMESON (1986) demonstrated possible utilization of this 
technology for geometrically complex problems. 

To summarize then, in the area of a posteriori error estimates, 
post-processing and adaptive mesh refinement, useful technology has been 
developed for linear structural analysis and certain nonlinear problems. 
However, further research and development including extensive numerical 
experiments will be needed to establish a methodology base for 
three-dimensional inelastic analyses involving complex loading histories. 

2.4 NONLINEAR ANALYSIS OF SOLIDS AND STRUCTURES 

Application of finite element methods to nonlinear structural analysis 
involves increasingly complicated geometrical and material models and pushes 
currently available supercomputing facilities to their limits. Codes 
specifically designed for supercomputers have demonstrated such computers' 
potential. Papers have appeared recently discussing performance improvement 
of finite element codes on vector and parallel machines utilizing algorithms 
tailored for a specific computer environment. BENSON, HALLQUIST (1986) 
discusses the rigid body algorithm implemented in the DYNA code which reduced 
the amount of computations significantly for explicit finite element time 
integrations. Implementation of an iterative solution algorithm based on the 
element-by-element preconditioner in the NIKE code is discussed by HUGHES, 
FERENCZ, HALLQUIST (1986). The utilization of supercomputers in an engineering 
environment involving the development of these codes is presented by GOUDREAU, 
BENSON, HALLQUIST, KAY, ROSINSKY, SACKETT (1986). 
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Investigations on vectori zabi 1 i ty of basic finite element operations continue. 
For instance, the details of vector-matrix multiplication in an 
element-by-element manner are discussed, including timing figures, for various 
finite element models by HAYES, DEVLOO (1986). An algorithm to take full 
advantage of vector machines in solving a system of ordinary differential 
equations is proposed by BROWN, HINDMARSH (1986). Iterative algorithms based 
on the ORTHOMIN method and incomplete factorization are studied by ZYVOLOSKI 
(1986). It is anticipated that the development of algorithms closely tied to 
the actual data manipulations and arithmetic operations in computing machines 
will eventually lead to faster finite element codes not only on supercomputers 
but also on smaller scalar machines. 

Utilization of parallel computing for equation solution is a subject 
investigated extensively in the last few years. A series of papers by Utku, 
Melosh and their collaborators looks into the possible parallelization of 
direct stiffness matrix factorization [UTKO, MELOSH, SALAMA, CHANG (1986); and 
UTKU, SALAMA, MELOSH (1986)]. The possible utilization of hypercube 
architecture is studied by NOUROMID, ORTIZ (1986) for the factorization and 
iterative solution of finite element equations. It is observed in these papers 
that a simple reduction in the total number of operations does not itself 
always provide the optimal algorithm for vector and parallel processing. An 
increasing number of sessions at scientific meetings are now being devoted to 
discussing this aspect of finite element computations, and extensive 
literature in this area is being accumulated. However, further research and 
development still needs to be done before robust algorithms become available 
for a wide range of engineering applications in order to utilize fully the 
potential of modern computing machinery. 

Notable progress has been made in the last few years in the field of 
computational plasticity, in particular, design and analysis of integration 
algorithms for rate-independent plasticity constitutive equations. ORTIZ, SIMO 
(1986) summarizes the development of such integration algorithms and 
investigates their accuracy in great detail. The basic idea is to generalize 
the return mapping algorithm based on the elastic predictor. Application to 
viscoplastic constitutive equations is also included in this paper. SIMO 
(1986) extends the results to finite deformation plasticity by developing a 
simple and efficient finite element algorithm in which elastic deformation of 
finite amplitude is accommodated. 

Approaches for capturing localization due to inelastic material response in 
global finite element computations have been developed which enable 
discontinuities to be modeled at the subelement scale. The work by HUGHES, 
SHAKIE (1986) extends the conventional return mapping algorithm for J2 flow 
theory to yield surfaces with corners. This extension enables the effects of 
localized plastic flow to be captured in a global manner. ORTIZ, LEROY (1986) 
proposes an algorithm which allows shear bands to generate inside a finite 
element. The numerical results indicate the potential of the method for 
capturing local failure of structures without excessive mesh refinement. 
Methods for bringing in local failure modes in finite elements to capture 
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localization effects under strain softening are discussed by WILLAM, SOBH, 
STURE (1986). A survey of the finite element technology for capturing 
localized failure modes is given by NEEDLEMAN (1986). This is a relatively new 
field and further work will be needed before localized microscopic material 
responses due to embedded singularities in large-scale structural systems can 
be calculated effectively. 

2.5 CONCLUDING REMARKS 

In the theoretical aspects of finite element development, the utilization of 
the mixed variational principle has been and continues to be the main thrust 
of research and development work. This approach not only provides tools to 
improve the performance of elements but also leads to a deeper insight on the 
solution algorithms for linear and nonlinear problems. The mixed iterative 
methodology developed under the HOST contract takes full advantage of these 
modern theoretical developments. Also, new ideas demonstrated in the MHOST 
program have attracted significant academic research interest in the 
computational mechanics community. 

In the application and computational aspects of finite elements, progress made 
in the last few years has not yet been fully incorporated into the MHOST code 
development. The information available in the literature indicates that 
further development to utilize modern algorithms and coding technology would 
enhance the performance of the MHOST code even more significantly. 
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3.0 GLOBAL SOLUTION STRATEGY AND ELEMENT FORMULATION 


3.1 INTRODUCTION 

This section summarizes the technology developed and implemented in the MHOST 
program during Task VB. 

The mixed variational formulation of Hu-Washizu type [HU (1955), and WASHIZU 
(1955, 1974)] is used as a driving mechanism to derive the iterative solution 
algorithms. A system of algebraic equations is generated, including the 
effects of the stress stiffening and the centrifugal mass terms in conjunction 
with the dynamic transient effects in a semi-discrete manner. The 
fully-discretized version of the finite element equation system is included in 
this report, with the Newmark family of algorithms used for temporal 
discretization. As demonstrated in ZIENKIEWICZ, LI, NAKAZAWA (1986), other 
time integration algorithms can be utilized to generate recursive forms for 
transient analyses. 

The implementation of iterative solution algorithms for the mixed finite 
element equations is briefly discussed in this section. A comparison of 
options made available in the MHOST code indicates that the secant version of 
the Davidon rank-one quasi-Newton update performs most reliably, except for a 
class of ill-conditioned problems. 

The updated Lagrangian formulation is used to facilitate large displacement 
effects. The inclusion of finite strain effects provides an experimental 
capability in the latest version of the MHOST code. A subsection is devoted to 
this development. 

Element formulations based on the selectively reduced integration and the 
assumed strain approaches are investigated. In the mixed iterative solution 
framework, element formulations applicable to anisotropic material response 
with rate independent deviatoric plasticity are identified as requiring 
further investigation. 

The subelement method is extended to handle inelastic response inside 
subelement regions. The validation cases included in this context indicate the 
performance of this new approach in capturing local inelastic responses around 
embedded singularities. 

3.2 SEMIDISCRETE FINITE ELEMENT EQUATIONS AND TEMPORAL DISCRETIZATION 

The finite element equations of dynamic equilibrium for a deformable body can 
be written in terms of the nodal acceleration vector a, the nodal velocity 
vector v and the nodal stress vector s as follows: 

M a + C v+ Bs = F, (1) 
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where Mis the mass matrix defined by 

M = P 

n 

with p and N being the material density and the vector of interpolation 
functions respectively. The damping matrix C may be defined as a sum of mass 
and stiffness matrices, M and K, such that 

C = PlM + P 2 K, (3) 

with P] and P£ being preassigned constants representing the effects of 
viscous and structural damping respectively. Symbol ical ly, the stiffness 
matrix can be written by 

K = / VN t D VN dx , (4) 

with V being the gradient operator and Dthe material modulus. The third term 
in Equation (1) represents the internal energy, with the B matrix being the 
discrete gradient operator defined in the finite element approximation' 
subspace such that 

B = 

12 

Here an equal order interpolation using the same basis functions is assumed 
for the acceleration, velocity, displacement, and stress. The nodally 
interpolated strain is recovered via 

e^G’ 1 B t u, (gj 

where G is the diagonalized gramrn matrix whose Ith diagonal entry is 
calculated by 

N p r 

G t = 2 / N t N, dx, 

1 J=1 1 J (7) 

with Np being the total number of nodes in a mesh. The nodal stress state is 
recovered by 


/ 


VN T N dx. 


(5) 




N t N dx 


( 2 ) 


s = D e . (3) 

Note that all the integrations indicated in the above equations are evaluated 
approximately using numerical quadrature. The particular choice of numerical 
integration rule is discussed in the subsection on element formulation. 
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Next, an abstract iterative process is constructed for the semi-discrete 
Equations (1) by means of displacement preconditioning. The approximation of 
internal energy by the product of stiffness matrix Kand total displacement 
vector u leads to the recursive expression 


M a-j+i +Cvi+] +KUj + ] s F- (B3, -KUj), ( 9 ) 

with subscripts denoting the iteration counter. For linear elastic problems, 
the displacement converges linearly given an appropriate initial condition 
such as u 0 = U. 

Note that the strain projection procedure and the stress recovery operation. 
Equations (6) and (8), are executed every time the displacement vector is 
updated. 

The fully discretized system of equations is derived for the mixed iterative 
solution by introducing the temporal approximation. For the sake of simplicity 
and clarity, the Newmark family of algorithms is introduced in a finite 
difference fashion, e.g., HUGHES (1984), such that dynamic equilibrium is 
satisfied at time t + At as follows: 

M a t+At + C v t+At + B s t+At = F t+At . ( 10 ) 

This then leads to a time discrete displacement preconditioning 


Ma* + f + Cv! + f + K u^ +At = F t+At - (Bs* +At - K u* +At ), 


i+1 


i+1 rrvu i+i 


with the updates for displacement and velocity given by 

u‘ +At = u‘ ♦ At v‘ + ^ {(!-») a‘ 4. 2P a‘ +At } , 


V t+ At = yt + At | (1 . 7) a t + 7 a^ +At | 


Nith the aid of 


Au^ +At = u^ +A ^ - u t 
1 1 


and 


gu t+ A t _ Au t+At - Au* +At , 


( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


a recursive form is obtained as follows: 
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where the incremental displacement Au^*^ is updated by 


Au 


t+At _ . t+At 
i+1 “ A i 


+ 


8u 


t+At 


(17) 


and the incremental strain is recovered at each node by 

Ae t+At _ q- 1 q T Au t+At . (18) 

In a general setting of history dependent inelastic constitutive integration, 
the stress recovery process is written as 

Ae t+At 

/ ’ sde 09) 

0 

which is fed into the equilibrium iteration defined by Equation (16). 

The computational effort in the above iterative process for direct time 
integration of the discrete equations of motion involves assembly of the 
element operator matrices which appear on the left-hand side of Equation (16). 
This requires exactly the same amount of effort as to solve the problem by the 
conventional displacement method. Additional computations required for the 
mixed solution involve a few back substitutions and strain and stress recovery 
at nodes, computations which are relatively inexpensive compared to matrix 
assembly and factorization. This additional computational cost would vanish 
when the scheme is applied to nonlinear problems in which matrix assembly and 
factorization or back substitution need to be performed repeatedly no matter 
which finite element method is used to drive the solution. Therefore, the 
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seemingly complicated derivation of the mixed iterative solution scheme does 
not increase the computational effort in comparison with the displacement 
formulation. Thus, the advantages of equal order interpolation of 
displacement, strain and stress can be fully utilized without fear of central 
processing unit (CPU) time penalties. 

With particular attention to turbo machinery blades and other rotating 
structures, the stiffness matrix can be modified to include the effects of the 
stress stiffening and the centrifugal mass terms. These results are based on 
linearized models of the large displacement and corresponding follower force 
equations under a given angular velocity. 

The centrifugal load vector is generated by integrating the body 

F Q s J N T po ; 2 r( x ) dx , 

Q 

witn w being the angular speed and r the vector distance between 
the axis of rotation. Under the centrifugal loading given above, 
iterative algorithm generates a certain stress field represented 
stress vector s 0 such that the discrete equilibrium equation 


is satisfied. Our interest here is to calculate the linearized response of the 
structure under an initial state of stress given by s 0 . If these initial 
stress terms are taken into account, the equilibrium equation system is 
modified (excluding the damping term) to 

M a + Kq u + B T B = F , (22) 

where K? is the geometric stiffness term associated with the initial stress 
field s 0 and is given by 


force 

( 20 ) 


a point and 
the mixed 
by the nodal 


K G = / V,SjT s ° VN dX * (23) 

Note that Equation (22) is a general expression for the linearized response of 
structural systems under a prescribed initial stress field. 

When this term, i.e.,KgU, is included in the mixed iterative computations, 
the abstract recursive form. Equation (9), becomes: 

M a 1 + l + (Kg + K ) u jtl * F - ( B s, - (Kg + K> Uj), (24) 
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or in terms of displacement update 


with 


(Kg + K ) Au i = F - ( B s i 


u 


i+1 


u . 


+ Au. 

l 


Kq u 1>* 


(25) 


In the residual force calculations, the initial stress matrix needs to be 
evaluated repeatedly with the stress field calculated at the beginning of the 
increment. 

The additional force due to the relative motion from the equilibrium position, 
referred to as the centrifugal mass term, is also added to the equilibrium 
equation system resulting in 

Ma + (Kg-M c )u+Bs = F , (26 

where M c is the centrifugal mass matrix defined by 

M c = J pJ N T L N dx , (27 

Q 


with L being a matrix based upon the components of the unit vector m parallel 
to the axis of rotation. L is given by 


L = 


2 2 

-m^m2 


~ m 3 m l 


2 2 
m^+mi 






2 2 
mi+m 2 


(28) 


When this term is activated i-n the quasi-static and transient dynamic 
computations, an additional correction of the residual vector calculation is 

made by replacing with (Kg - M c ). 

Typically, these optional terms are activated for the vibration analysis of 
rotating structures by extracting the modes of the eigenvalue problem 

|m-\(Kq-M c +K)|x-0. 
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In the modal analysts, the displacement preconditioner is assumed to represent 
the structural stiffness. However, the stress values used in the evaluation of 
initial stress terms are calculated by the mixed iterative method. Thus, a 
more accurate stress field is generated than via the conventional displacement 
method. Hence, the solution of Equation (29) is expected to be somewhat more 
accurate than that associated with the displacement method approach. 

For details of derivations of initial stress and centrifugal mass matrices and 
fully worked out examples, see THOMAS, MOTA SOARES (1973) and RAWTANI, 
DOKAINISH (1971). 

In summary then, Equation (16) represents the recursive form of the basic time 
discrete displacement preconditioning as given by Equation (11); furthermore, 
Equation (11) itself is a discretized form of the original basic equation of 
dynamic equilibrium as given by Equation (1). The bottom line, however, is 
that the sequence of computations for the basic equation set is given by the 
order of Equations (16) through (19). For centrifugal load effects, Equations 
(24) and (25) replace Equations (16) and (17) and Equations (20) and (23) are 
also added. 

3.3 GLOBAL SOLUTION ALGORITHMS 

In the present implementation of the mixed iterative solution procedures, the 
displacement stiffness equations need to be assembled and factorized. To 
enhance the efficiency of the MHOST code, a new profile solver algorithm has 
been developed to replace the constant bandwidth Crout decomposition 
algorithm. The implementation of this new profile solver is based on code 
published by TAYLOR (1985). Examples presented in that paper include both 
symmetric and nonsymmetric cases using Crout decomposition. Implementation of 
the profile solver into the MHOST code excludes a nonsymmetric capability in 
order to gain maximum efficiency. The profile storage of the global stiffness 
equations reduces memory requirements as well as time required for 
factorization. Additional computational overhead is required to prepare the 
elimination table for a stiffness matrix stored in profile form. A table 
comparing the performance of these solvers is included in Section 4.4.1. 

In the quasi -static computations, a series of options to accelerate the rate 
of convergence for nonlinear iterations has been made operational with the new 
profile solver capability. These options are as follows: the modified Newton 
method, the quasi-Newton method of inverse BFGS update, the secant Newton 
implementation of Davidon rank-one quasi -Newton update, the conjugate gradient 
method, and the line search method. The spherical path version of the arc 
length method incorporating modified Newton iteration is also now available 
with the profile solver capability. For details of the implementation of these 
numerical processes in the mixed iterative solution framework, see the third 
year progress report [NAKAZAWA (1986)]. 

The eigenvalue extraction procedure implemented into the MHOST code is based 
on the subspace iteration approach documented by BATHE, WILSON (1976). In this 
process, the large eigenproblem of structural systems is mapped into a 
subspace of finite dimensions and Jacobi iteration is performed to solve the 
small eigenproblem set in the subspaces. This eigenvalue extraction procedure 
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is used to compute the natural frequency and vibration mode of unstressed and 
prestressed structures as well as the buckling analysis of prestressed 
structures undergoing elastic or inelastic deformation of infinitesimal or 
finite amplitude. 

Note here that in the eigenvalue extraction, the standard displacement 
stiffness matrix is used to represent the response of structural systems. 
However, in buckling and modal analyses with prestress, the improvement of the 
stress field generally improves the quality of the solution considerably. The 
utilization of mixed stress interpolation in the eigenvalue analysis can be 
viewed as equivalent to a perturbed eigenvalue problem in which the effect of 
the independent stress approximation is regarded as the perturbation to the 
displacement finite element system. It is anticipated that an algorithm 
proposed by SIMO (1985) and implemented for the probabilistic finite element 
code at MARC may be usable for eigenvalue extraction of mixed systems of 
finite element equations. 

3.4 LARGE DISPLACEMENT ANALYSIS 

In this section, the large deformation algorithm implemented in the mixed 
iterative method is described. The formulation utilizes an updated Lagrangian 
mesh approach. The choice of deformation and stress measures in finite 
deformation analysis depends primarily on the class of material involved, 
along with the necessity that the measures are objective (invariant) for rigid 
body motions. For elastic-plastic behavior where the response depends 
primarily on the current state of stress, a Cauchy stress and 
rate-of-deformation (i.e., strain rate) constitutive formulation is frequently 
most convenient. 

Beyond the choice of a Lagrangian mesh description, the formulation can be 
further simplified by evaluating the equations of motion at the current 
(deformed) configuration. This is achieved by continually updating the mesh 
geometry to the current state. Utilizing an updated geometry results in 
simplification because all matrix expressions take the same form as in small 
deformation theory. The only additional quantities needed are deformation 
gradients and rotation tensors evaluated at nodes, along with follower force 
matrices. This formulation is referred to as the updated Lagrangian approach 
in the literature. 

The updated Lagrangian algorithm for large deformation problems is essentially 
the same as that for nonlinear small deformation problems, with the exception 
that all tensors and integrals are evaluated with respect to the current 
configuration. The updated Lagrangian algorithm presented in this section 
involves nodal coordinates being updated continuously and deformation 
gradients being evaluated at nodes. 

Table 1 presents a summary of the algorithm for an increment. All new terms in 
this table are defined in the Nomenclature section. Following initialization, 
the process is repeated iteratively until equilibrium has been satisfied. 
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Table 1 Updated Lagrangian Algorithm 


Initialization : 

u 1 = u 1 * 1 [a] 1 ' = [a] 1 '" 1 R 1 = F 1 

Equation Formulation and Displacement Solution : 

Au 1 = ( K 1 )" 1 R 1 


Update Geometry : 


x i+1 = x 1 + s 1 Au’ 

pi+1/2 

rel 

= I + s 1 v(| Au 1 ) 

u i+1 = u 1 + s 1 ' Au 1 


= I + s 1 V( Au^ ) 

pi+1/2 _ pi+1/2 pi 

rel 

R it,/2 

= pi+1/2 (u i+l/2 } -l 

F i + 1 = F 1+J F i 
rel 

R i + 1 

= F 1+1 ( U i+1 ) _1 

Strain Projection: 



Large Strain 

A L ej i+1 = (CS 1 + 1 f T B 1 + 1 s 1 Au 1 

Small Strain 

A[e] 1+1 = ( CSV 1 B 1 s 1 Au 1 


Stress Recovery : 

[<xjj +1 = D 1 + 1 ( ROT 1 + 1/2 ) T A[€] 1 + 1 ROT 1 + 1/2 + [o’] 1 
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Table 1 Updated Lagrangian Algorithm (continued) 


Rotate Back to Global System : 

= ROT 1 + 1 |>] | J +1 ( ROT 1+1 ) T [e ]^ 1 = ROT 1 + 1 [e]j+ 1 (ROT 1+1 ) T 

R 

D 1 + 1 = ROT 1+1 D 1 ^ 1 (ROT i+1 ) T [o !] i+1 = ROT i+ 1 [a]’ +, ( ROT i+1 ) T 

Form Residual : 

R i+1 =/ (N i+1 ) T f i+1 N i+1 dx 1 + 1 

+ / ( N 1+1 ) T h 1+1 N 1+1 dS 1+1 

an 1+1 

- (B i+1 ) T |>] i+1 

If converged, exit; otherwise update displacements via line search or 
displacement solution and repeat. 
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In the finite displacement calculations, the stiffness array is assembled in 
the following manner and repeatedly updated in the iteration loop unless 
another iteration process (modif ied-Newton, quasi-Newton, or secant-Newton 
option) is specified. Note that the stiffness matrix does not correspond to 
the classical Newton method in the mixed iterative solution framework because 
of the independent C°-continuous stress interpolation functions. Here, B is 
the conventional strain displacement matrix. 


K 1 = / B T ( x 1 ) D| B(x 7 ) dx 1 

n 1 


Tangent stiffness matrix 


+ / B^x’JO^Kx 1 ) dx 1 

Initial stress matrix 


- / n (x 1 ) f 1+ Nix 1 ) j] I J 1 )" d| 


- / N^x 1 ) h 1+1 Nix 1 ) J 1 , ( J 1 ) -1 d£ 

8* * 


Body force follower force 
due to volume change 

( 30 ) 


Surface traction follower 
force due to surface area 
change 


- f [vnV) f^Ntx 1 ) + N T ( x 1 ) 

oi 


f 1 + 1 VN(x 1 ) ] dx 1 

Body force follower force 
due to rotation 



[VN^x 1 ") h 1 + 1 N(x 1 ) + nV) h 1 + 1 VN(x 7 )] dS 


Surface traction follower 
force due to rotation 


Whenever the stiffness matrix is reformulated in an updated Lagrangian 
algorithm, the shape functions, derivatives of the shape functions and 
determinants of Jacobians are evaluated at the nodal coordinates of the 
current configuration. In addition, the material and stress matrices used are 
those computed in the immediately preceding iteration. 
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4.0 COMPUTER PROGRAM DEVELOPMENT AND VALIDATION 


4.1 INTRODUCTION 

Version 4 of the MHOST code consists of over 47,000 lines of FORTRAN source 
statements including extensive comments. The concept of multiple-analysis 
drivers has been used to maintain the quality and efficiency of the code. A 
brief discussion of the architecture of the code is presented in this report. 

Note that the multiple driver strategy has kept the code usable for well 
tested capabilities, while new options and enhancements of existing features 
were being developed and tested. Also, such a strategy has helped to keep the 
code readable since the number of conditional statements is minimized. In 
order to use this code effectively for production purposes in engineering 
environments, further streamlining of the code would be needed to improve the 
performance of this finite element package. 

The user interface has been enhanced to deal with numerous algorithmic options 
added to Version 4 of this code. A subsection briefly summarizes the program 
features and user commands. 

All the development work has been carried out on a PRIME 9955 at MARC Analysis 
Research Corporation. The PRIME FORTRAN 77 compiler is employed without 
extensive utilization of extensions to ANSI standard. Versions for VAX/VMS, 
IBM/CMS and CRAY/COS environments have been created and tested using their 
FORTRAN 77 compilers. The utilization of FORTRAN 77 features which are not 
supported by the previous FORTRAN standard have made the new code incompatible 
with the old FORTRAN compilers. 

A profile solver approach has replaced the constant bandwidth Crout 
decomposition routines. The new solution package improves the efficiency of 
the code significantly in terms of both core storage and computing speed. A 
large portion of the MHOST code has been rewritten to take advantage of this 
new solver capability, including the quasi-static analysis, eigenvalue 
extraction and transient time integration subsystems. The out-of-core frontal 
solution subsystem remains unchanged and may be employed for large-scale small 
deformation quasi-static calculations. 

4.2 COMPUTER CODE ARCHITECTURE 

The architecture of the MHOST code is schematically shown in Figure 1. The 
execution supervisor routine (SUBROUTINE HOST) controls a multiple number of 
analysis modules in a consistent manner. In Version 4 of the MHOST code, a 
pair of subroutines has been added to check the consistency of the Parameter 
Data and to generate internal flags for the selection of analysis modules. The 
objective here is to avoid combining features of the program package not 
intended by the developers to be combined. 
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The analysis modules provide the control structure for the incremental 
iterative algorithms implemented in the MHOST code. The source program with 
comments is designed to serve as the schematic flow chart of the computational 
process. 

The schematic flow of the element and nodal data manipulations is coded in the 
element assembly submodules which are entered from the analysis modules. The 
element assembly submodules perform operations independent of element type and 
constitutive model. The assumptions introduced in the mechanics aspects of the 
formulations are explicitly coded at this level. These submodules call the 
librarian routines for elements and constitutive models. 

The element librarian subprogram (SUBROUTINE OERIV) generates quantities 
unique to the element used in an analysis. The MHOST code uses for the most 
part the standard finite element matrix notation, as in ZIENKIEWICZ (1977). 

The element specific information returned from, the librarian subprogram is the 
strain-displacement array referred to as the B matrix in previous subsections. 

The constitutive equation librarian subprogram (SUBROUTINE STRESS and BMSTRS) 
is designed to accommodate the nodal storage of stress and strain values. The 
STRESS subroutine sets up the loop over the points at which constitutive 
equations are evaluated. The current implementation is a nested double loop 
with the outer loop being over the nodes and the inner loop over the 
integration layers through the thickness. Note that the conventional 
displacement method can be recovered by restructuring this loop in conjunction 
with a few minor modifications in the core allocation for stresses and strains. 

The librarian subprogram calls the constitutive equation package (SUBROUTINE 
NODSTR) from which individual subprograms for initial strains, stress recovery 
and material tangent are accessed. The librarian subprogram also controls the 
pre- integration of stresses, strains and material tangent over thickness for 
the shell element. 

The evaluation of the constitutive equations is one of the most costly 
operations in the nonlinear finite element computations. An attempt has been 
made to minimize the execution of this process during the incremental 
iterative analysis, typically once every iteration during the recovery of the 
residual vector. At the beginning of each increment, this process may be 
executed to evaluate the material tangent as well as the contribution of 
initial strain terms which are often necessary for proper displacement 
pre-conditioning. 

Except for a small amount of information related to the convergence of the 
iterative solution, all report generation is performed at the end of an 
increment. Optionally, post-processing and restart files are written at the 
end of user-specified increments. Generic reporting subprograms are called 
from analysis modules inside the loop over the increments. 
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A brief description of major subprograms is provided below. 
Execution Supervisor 


MAIN PROGRAM - declares the work space in blank common as an integer 

array. Also defines system parameters which are machine 
independent. Then passes control to the actual execution 
supervisor SUBROUTINE HOST. 

SUBROUTINE HOST - controls the sequence of execution of analysis modules. 

First, this routine executes the control parameter data 
reader SUBROUTINE DATIN1 and checks for consistency by 
entering SUBROUTINE LETCMD and RUNCMD. The current 
structure of the code allows certain combinations of two 
analysis modules to be executed sequentially. Analysis 
modules called from the execution supervisor are as 
follows: 

SUBROUTINE STATIC for a quasi-static incremental iterative 
solution. 

SUBROUTINE DYNAMT for a transient time integration of the 
dynamic equilibrium equation system in an incremental - 
iterative manner. 

SUBROUTINE MODAL for eigenvalue extraction in vibration 
mode analysis. This subsystem may be executed after a 
quasi-static analysis for modal analysis of prestressed 
structures. 

SUBROUTINE BUCKLE for eigenvalue extraction in buckling 
load calculations. This subsystem is executable only after 
a quasi-static analysis. 

SUBROUTINE FRONTS for a quasi-static incremental iterative 
solution by the out-of-core frontal solution procedure. 
Note that certain options are not available in this 
subsystem. 

SUBROUTINE SUPER for a linear dynamic response calculation 
by the method of mode superposition. 


Input Data Reader 

There are three major subprograms: 

SUBROUTINE DATINl - Reads and interprets the parameter data input called by 

the execution supervisor. All the default values for 
control variables are set in this routine. 


22 



SUBROUTINE BULKIN - Reads and interprets the finite element model definition 

data and prints the mesh and loading data for the initial 
increment (number 0). The bulk data reader is entered 
before the execution of analysis modules. This subroutine 
utilizes the following lower level routines: 

SUBROUTINE INITI1 for the memory allocation of integer 
work space in the blank common to store nodal and element 
data. 

SUBROUTINE 0ATIN2 for the model data input. 

SUBROUTINE DAT0U1 for reporting the model definition data. 

SUBROUTINE CHKELM for the detection of clockwise element 
connectivity which results in a negative Jacobian matrix. 
This routine automatically corrects the connectivity table 
to a counterclockwise direction. 

SUBROUTINE SUBDIV for memory allocation of subelement mesh 
data and automatic mesh generation of subelements. 

SUBROUTINE INCRIN - Reads and interprets the loading and constraint data for 

each increment. This routine is invoked by individual 
analysis modules from inside the loop over the increments. 

SUBROUTINE DATIN3 - A small subset of the bulk data reader SUBROUTINE DAT INI , 

is used for actual operations including initialization of 
arrays at the beginning of an increment. 

Algebraic Operations Subsystem 

There are three packages of routines for the profile solver, frontal solver 
and eigenvalue extraction. The profile solution package consists of: 

SUBROUTINE COMPRO - Sets up the integer array for the profile of the global 

stiffness equations to be stored in a profile form. 

SUBROUTINE ASSEM5 - Assembles the element stiffness equations into the global 

equation system stored in a profile form. 

SUBROUTINE S0LUT1 - Controls the iterative solution processes including the 

vector update required for the quasi- and secant-Newton 
iterations. 

SUBROUTINE DECOMP - Factorizes the global stiffness equations stored into a 

profile form. 

SUBROUTINE SOLVER - Performs the back substitution and generates the update 

vector for the incremental displacement. 
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The frontal solution package consists of: 

SUBROUTINE FRONTW - Estimates the front matrix size to be accommodated in the 

core memory. 

SUBROUTINE INTFR - Allocates memory for the work space required for the 

frontal solution. 

SUBROUTINE PRFRNT - Sets up the elimination table for the frontal solution. 

SUBROUTINE FRONTF - Assembles and factorizes the global stiffness equations 

simultaneously. 

SUBROUTINE FRONTB - Performs the back substitution and generates the updates 

for the displacement vector. 

SUBROUTINE FRONTR - Controls the iterative solution processes including the 

vector update required for the quasi- and secant-Newton 
iterations. Calls SUBROUTINE FRONTB for the incremental 
displacement update vector. 

SUBROUTINE VDSKIO - Controls the data stream stored in the in-core buffer area 

and the actual out-of-core storage devices. 

The eigenvalue analysis package consists of: 

SUBROUTINE EIGENV - Controls execution of the eigenvalue extraction subsystem. 

SUBROUTINE INIMOP - Initializes the array for eigenvectors. 

SUBROUTINE SUBSPC - Performs the subspace iteration and generates a specified 

number of eigenvalues and eigenvectors. 

SUBROUTINE JACOBI - Solves the eigenvalue problem in the subspace by Jacobi 

iteration. 

There are a number of subprograms used commonly by the Algebraic Operations 

Subsystem: 

SUBROUTINE STRUCT - Controls memory allocation for global algebraic 

manipulations at the beginning of every increment. 

SUBROUTINE INITI2 - Allocates memory required for storage of the global 

stiffness matrix and other vectors necessary in the linear 
algebraic manipulation of the finite element equations. 

SUBROUTINE LINESR - Calculates the search distance when the line search option 

is turned on. 
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Element Assembly Submodules 


These are subprograms constructing vectors and matrices appearing in the 
algorithmic description of the mixed iterative process discussed in previous 
subsections. 


SUBROUTINE ASSEM1 - Assembles the displacement stiffness matrix for 

preconditioning purposes. All the options for kinematic 
and constitutive assumptions are tested in this module. 

SUBROUTINE ASSEM2 - Assembles the coefficient matrix for the transient time 

integration by the Newmark family of algorithms. This 
routine is evolved from SUBROUTINE ASSEM1 and contains all 
the options. 

SUBROUTINE ASSEM3 - Assembles the coefficient matrix for the quasi-static 

analysis using the frontal solution subsystem. Large 
displacement, stress stiffening and centrifugal mass terms 
are not available in this package. 


SUBROUTINE ASSEM4 - Calculates the nodal strain and recovers the residual 

vector in a mixed form. The subelement solution package is 
entered from this subprogram. 


Element Loop Structure and Library Routines 


In the element assembly submodules, element arrays are generated in loops over 
the elements. The protocol for accessing the element library involves a 
sequence of subroutine calls as follows: 

SUBROUTINE ELVULV - Sets up the current element parameters (see Table 2 for 

variables and values) from the element library table. 

SUBROUTINE CNODEL - Pulls out quantities for the current element from the 

global nodal array and restores them in the element work 
space. Coordinate transformations necessary for beam and 
shell elements are performed in this subprogram. 

SUBROUTINE DERIV - Sets up the displacement-strain matrix for the current 

element by calling the element library subroutines. Those 
are: 

SUBROUTINE BPSTRS for plane stress elements, types 3 and 

101 . 


SUBROUTINE BPSTRN for plane strain elements, types 11 and 

102 . 
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SUBROUTINE BSOLID for three-dimensional solid element, 
type 7. 

SUBROUTINE BSHELL for three-dimensional shell element, 
type 75. 

SUBROUTINE BAXSYM for axi symmetric sol id-of -revolution 
elements, types 10 and 103. 

SUBROUTINE BTBEAM for linear Timoshenko beam element, type 
98. 

SUBROUTINE BASPST for the assumed stress plane stress 
element, type 151. 

SUBROUTINE BASPSN for the assumed stress plane strain 
element, type 152. 

SUBROUTINE BASSOL for the assumed stress three-dimensional 
solid element, type 154. 

SUBROUTINE UDERIV - A slot for a user coded element B matrix routine. 

The following subprograms are used to calculate terms appearing in the finite 
element equations: 

SUBROUTINE LMPMAS - Calculates nodal weight factors for the strain projection. 

SUBROUTINE STIFF - Performs matrix triple products to assemble the element 

stiffness matrix and the element load vector associated 
with the initial strain terms. 

SUBROUTINE STRAIN - Calculates element strains at specified sampling points 

and projects to nodes. 

SUBROUTINE CNSMAS - Assembles the consistent mass matrix for modal and 

transient analysis. 

SUBROUTINE INITST - Generates initial stress terms for quasi -static, buckling 

and modal analysis of prestressed structures. 

SUBROUTINE CENMAS - Evaluates the centrifugal mass terms for rotating 

structures at speed. 

SUBROUTINE RESID - Calculates the element residual vector for the global 

element. 

SUBROUTINE SUBFEM - Performs the subelement solution and calculates the 

element residual vector for the subelement mesh. 
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SUBROUTINE RELDFG - Calculates the relative deformation gradient at the 

element sampling point and projects to a node. 

SUBROUTINE RESDYN - Calculates the contribution of mass and damping terms in 

the element residual vector when the transient dynamics 
option is used. 

Material Library 

A system of subroutines is included in the MHOST code which covers a wide 

range of material models and initial strain assumptions: 

SUBROUTINE SIMPLE - Integrates the stress over the increment assuming the 

elastic-plastic response of the material is represented by 
a total nonlinear elastic secant modulus. Also generates 
the material modulus matrix. 

SUBROUTINE PLASTS - Integrates the stress over the increment and calculates 

the plastic strain by using the radial return algorithm. 

SUBROUTINE PLASTD - Calculates a consistent elastic-plastic modulus if the 

incremental equivalent plastic strain is positive. 

SUBROUTINE WALKEQ - Integrates the Walker unified creep plasticity 

constitutive equation and also generates a material 
modulus based on the temperature dependent elasticity 
assumption. 

SUBROUTINE LELAST - Calculates the stress for the constant material modulus 

given as data. 

SUBROUTINE THRSTN - Calculates the thermal strain. 

SUBROUTINE CRPSTN - Calculates the creep strain. 


4.3 USER INTERFACE 

A large number of options are made available in the MHOST code. A user selects 
options by specifying certain keywords in the parameter data section of the 
input data file. Some of the infrequently used features available in the 
previous versions are deleted for the sake of clarity. Table 2 summarizes the 
analysis capability of MHOST Version 4.2. 


i 

i 
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Table 2 MHOST Analysis Capability 


Element 




Three- 

Three- 

Definition 

Plane 

Plane 

Ax i symmetric 

Dimensional 

Dimensional 

Options Beam 

Stress 

Strain 

Solid 

Solid 

Shell 

Li near 






Isotropic 
Elasticity x 

X 

X 

X 

X 

X 

Ani sotropi c*! 
Elasticity 

X 

X 

X 

X 


Composite*^ 
lami nate 





X 

Simplified 

Plasticity 

X 

X 

X 

X 

X 

Elasto-*^ 
PI astici t y 

X 

X 

X 

X 

X 

Uni fied 

Creep-Plasticity 

X 

X 

X 

X 

X 

Thermal *2 
Strain 

X 

X 

X 

X 

X 

Creep*^ 
Strai n 

X 

X 

X 

X 

X 

Large 

Displacement x 

X 

X 

X 

X 

X 

Finite 

Strain 

X 

X 

X 

X 



Notes : 
*1 
*2 


Applicable only to linear elasticity. 

Not applicable to the unified creep-plasticity model in which the 
quantities are the integrated part of the model. 


*3 Anisotropic yield surface can be defined for a continuum element. 
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Table 2 MHOST Analysis Capability (continued) 


Analysis 

Module 

Option 

Beam 

Plane 

Stress 

Plane 

Strain 

Axi symmetri c 
Solid 

Three- 

Dimensional 

Solid 

Three- 

Dimensional 

Shell 

Quasi -static 
Analysis 

X 

X 

X 

X 

X 

X 

Buckl ing 
Analysis 

X 

X 

X 

X 

X 

X 

Modal 

Analysis 

X 

X 

X 

X 

X 

X 

Transient 

Dynamics 

X 

X 

X 

X 

X 

X 


In the model data section of the input data file, certain nonstandard features 
are made available. This subsection summarizes program features available in 
Version 4.2 of the MHOST code by the keywords. For full details, see the MHOST 
User's Manual . 

*AN ISOTROPY 

The anisotropic material property option is flagged by this parameter. The 
orientation vector for the material axes must be added in the model data 
section. The material properties along the material axes must be given by 
either the *DMATRIX option (for continuum elements, type 3, 7, 10 and 11 ) or 
the ^LAMINATE option (for shell element type 75). The anisotropic plastic 
response of a material is described by the user subroutine ANPLAS as 
documented in the MHOST User's Manual. 

*BFGS 

The inverse BFGS update procedure is invoked by flagging this option parameter 
with the default iterative algorithm being the straightforward Newton-Raphson 
scheme. 

♦BOUNDARY 

The nodal displacement constraints are imposed by virtue of penalization when 
the frontal solution subsystem is invoked. 
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*BUCKLE 


The initial stress matrix is generated and the eigenvalue extraction is 
performed to obtain buckling modes. This option can be invoked at an arbitrary 
step of the incremental non linear solution process so as to detect the change 
in the buckling load due to the inelastic response of the structure. 

♦COMPOSITE 

This parameter data line invokes the composite laminate analysis option for 
shell element type 75. Added to the MHOST code upon request from the Contract 
Monitor at NASA-Lewis Research Center. This option resets the element 
parameters. The number of integration layers is reduced to 1 and all the 
components of stress resultant are used as nodal variables. In the model data 
section, the user has to supply all the components of the stress-strain matrix 
by invoking the laminate option. 

♦CONJUGATE-GRADIENT 

The conjugate-gradient method is invoked by adding this parameter. The line 
search option is automatically turned on. Note that this option cannot be used 
with other iterative algorithms such as the BFGS or secant-Newton methods. 

♦CONSTITUTIVE 

Three constitutive approaches are available for describing material behavior, 
including the secant elasticity procedure (simplified plasticity) in which the 
material tangent is generated for Newton-Raphson type iterative algorithms, 
the von Mises plasticity procedure with the associated flow rule treated by 
the radial return algorithm, and the Walker nonlinear viscoplastic model in 
which an initial stress iteration using elastic stiffness is utilized. For 
experimental purposes, the linear elasticity option can be flagged, with the 
default option being the conventional plasticity model. 

♦CREEP 

The creep effect is taken into account via the time history being integrated 
in an explicit manner under control of an optional self adaptive time-step 
size algorithm. 

♦DISPLACEMENTMETHOD 

This option invokes the conventional displacement model in which the residual 
loads are evaluated directly at integration points. For linear elastic stress 
analyses, no iteration would take place when this parameter card is flagged. 

In inelastic analyses, the material tangent is interpolated and multiplied by 
integration point strains directly at each quadrature point. This option 
cannot be used for the creep-plasticity model which does not generate correct 
material tangent information at integration points. 
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*DISTRIBUTELOAD 


The body force and surface traction loadings are referred to as the 
distributed loads in the MHOST program. The body force option includes gravity 
acceleration definable in any direction, as well as centrifugal loading with 
the centerline and angular velocity being user-specified. 

♦DUPLICATENODE 

The nodal points for stress recovery can be disconnected by defining two nodal 
points at the same geometrical location and connecting them by this option 
which assumes the compatibility of displacement at these points. This option 
is used to define generic modeling regions and their interconnections. 

♦DYNAMIC 

The generalized Newmark solution algorithm is entered by using this parameter 
card. A primitive version of the adaptive time-stepping algorithm is 
incorporated if requested by the user. 

♦ELEMENTS 

The elements available in this version of the code are summarized in Table 3. 
Core allocation is performed for both nodal- and element-quantities based on 
the maximum storage space requirements among the types of element specified in 
this option. All the element types must be specified here, including those 
only associated with the subelement regions. 

♦EMBED 

The subelement iteration technology is flagged by this option to signal the 
code to allocate working storage for subelement data in a hierarchical manner. 

The actual subelement mesh definition and the nodal- and element-data storage 
allocation take place when the individual subelement is defined. 

*FINITE_STRAIN 

This keyword invokes the finite-strain option at a specified load/time 
increment. 

*FOLLOWER_FORCE 

The loading data is assumed as the follower force when this keyword is 
included in the parameter data section. 

♦FORCES 

Concentrated nodal forces are defined and stored in an incremental manner. The 
core allocation takes place only when this option is invoked. 
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*FRONTALSOLUTION 


The frontal solution option for quasi -static analysis is available in this 
version of the code and utilizes out-of-core storage devices, thus increasing 
the capacity of the program significantly. 

*HARDENING_SLOPE 

The maximum number of break points for a piecewise linear approximation of 
uniaxial elastic-plastic behavior is specified with this parameter keyword. 

*LARGE_DI$PLACEMENT 

This keyword invokes the large displacement option at the specified load/time 
increment. 

*LINE-SEARCH 

The line search option requires positive action by the user to be turned on. 
This algorithm is usable in conjunction with all iteration algorithms 
available in the MHOST program. When the conjugate-gradient iteration is used, 
the program automatically turns on the line search option. 

*LOUBIGNAC 

Parameters for numerical quadrature used in mixed iterative processes are 
definable in a very precise way. To construct the stiffness matrix the full, 
selective, or selective with filtering scheme can be chosen. For the residual 
vector integration, full or reduced integration can be chosen. The strain 
integration can be performed by using either uniformly reduced integration, 
trapezoidal integration with a reduced shear strain approximation, or 
trapezoidal integration with a filtering option. 

*M0DAL 

Free vibration modes are extracted by invoking this option for linear elastic 
structures only. The subspace iteration technique is utilized together with 
the power shift option. 

*N0DES 

All the variables are defined and reported at nodal points. In the incremental 
processes, the deformation and stress history are stored only at the nodal 
points. Note that this architecture economizes storage substantially as 
compared to the fully integrated finite element displacement methods. 

*N0ECH0 

This option suppresses the echo print of the model data. 
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♦OPTIMIZE 


Bandwidth optimization based on the CuthiTl -McGee algorithm is available in 
the MHOST code, whereas no optimization for the frontal solver is available. 

♦PERIODICLOADING 

For transient analyses, nodal displacements and concentrated forces can be 
given as sinusoidal functions via this option; this reduces the effort needed 
to prescribe the load history for a class of simple test problems. 

♦POST 

The post-processing file, which contains all the information supplied to and 
produced by the code, is generated on the file connected to FORTRAN unit 
number 19. This file is formatted and can easily be manipulated by 
commercially available post-processing packages with minor modifications. A 
version of MENTAT supports the MHOST post-processing file. The header record 
of the file is designed to be compatible with the MARC post file, details of 
which are published and processed by many finite element graphics packages. 

♦PRINTSETS 

Report generation is carried out on a nodal point basis, with the element 
integration point option supported by interpolation processes using 
appropriate shape functions. 

♦REPORT 

The frequency of the line printer report generation is now controlled by this 
option with the default being to print at every increment. 

♦RESTART 

This parameter card invokes the program to read the restart tape and to set up 
the system to resume the analysis from the point where the tape is written. 

♦SCHEME 

Parameters for defining characteristics of the time integration operator can 
be specified by the user. The default is the average acceleration algorithm 
commonly used in nonlinear dynamic finite element analyses. 

♦SECANT-NEWTON 

This parameter activates the secant-Newton implementation of the Davidon 
rank-one quasi -Newton algorithm. The core-storage requirements of this 
procedure are significantly less than for the BFGS update and the execution is 
relatively fast. This parameter is recommended for inelastic analyses of 
sol id-continua. 
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♦STIFFENING 


This keyword is used to invoke initial stress terms for quasi -static, 
transient dynamic and modal analyses. 

♦STRESS 

As an option, the boundary condition for stress can be specified by the user, 
although no mathematical justification is yet available for this type of 
boundary loading. Any stress component can be prescribed at any nodal point. 
Simple numerical tests have shown that inconsistent imposition of boundary 
stress values indeed lead to rapid divergence in the iterative process. 

♦TANGENT 

In the context of modified Newton iterations, the tangent stiffness matrix for 
subsequent iterations is user-definable. When this option is invoked, the 
default is the first tangent matrix used during the current increment, not the 
elastic stiffness matrix. This type of implementation is referred to as the 
Kyi method of modified Newton iteration. In the BFGS process, the same 
approach is used and no options are made available. In the modified Newton 
iteration, the tangent array is fixed after a specified number of updates. 

♦TEMPERATURE 

Nodal temperatures are read as input and used to generate thermal strains. 
These temperatures are used also during the evaluation of creep strains and 
the integration of the coupled creep plasticity model (Walker's model). 

♦THERMAL 

Temperature dependent material properties are evaluated when this option is 
invoked, and an appropriate user subroutine has to be provided to the system 
prior to execution. This operation is not necessary for the creep plasticity 
model in which the temperature dependency is readily incorporated. 

♦TRANSFORMATIONS 

Coordinate transformations at nodal points are specified by this option in 
which the angle of rotation is provided by the user so that the code can 
generate necessary transformation matrices. The post-processing file does not 
support the coordinate transformations. 

♦TYING 

Multiple degree-of -freedom constraint equations are specified by the user 
through this option. Under certain circumstances, this option is slightly more 
flexible than the duplicate node option; for instance, the user can disconnect 
three elements at a point. Note that the constraint equations are generated 
only for displacement degrees -of-freedom, which are calculated only in an 
implicit manner. 
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The following parameters are used to signal to the MHOST program the presence 
of specific user subroutines: 

*UB0UN 

*UC0EF 

*UDERIV 

*UF0RCE 

*UH00K 

*UPRESS 

*UTEMP 

*UTHERM 

*WK$LP 

As summarized above, the analysis capabilities included in the MHOST code 
cover most requirements for calculating inelastic deformations of turbine 
engine hot section components. The MHOST code has evolved into a General 
Purpose Structural Analysis package with strong emphasis on nonlinear material 
behavior usable for a wide range of problems. The free format data input 
routines and report generation packages have been improved in the Task II 
program development effort so as to create a comfortable environment for code 
users. 

Table 3 summarizes the elements currently available in the MHOST Version 4.2 
code and lists parameters internally used to define element characteristics. 
Further detail is available in the MHOST User's Manual. 
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Table 3 Element Parameters 


Beam P. Stress P. Strain Axsym. Brick Shell 


ITYPE 

NELCRD 

NELNFR 

NELNOD 

NELSTR 

NELCHR 

NELINT 

NELLV 

NELLAY 

NDI 

NKHEAR 

JLAW 


NELCRD Number of coordinate point data values per node. 

NELNFR Number of degrees -of-freedom per node. 

NELNOD Number of nodes per element. 

NELSTR Number of stress and strain components per node. 

NELCHR Number of material property data values for the element. 

NELINT Number of 'full' integration points per element. 

NELLV Number of distributed load types per element. 

NELLAY Number of layers of integration through the thickness of the 
shell element. 

NDI Number of direct stress components. 

NSHEAR Number of shear stress components. 

JLAW Type of the constitutive equation. 
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4.4 VALIDATION 


4.4.1 Performance Improvement by the Profile Solver 

Table 4 summarizes computing time and memory requirement reductions after 
implementation of the profile solver solution package into the MHOST code. The 
performance figures are applicable for the MARC Corp. PRIME 9955 computer 
using a Primos FORTRAN 77 compiler with the full optimization option turned 
on. The number of words appearing in the table apply to situations involving 
single precision (32 bits) words in the blank common work space. For 64 
bits /word machines such as CRAY, the memory requirement in terms of number of 
words will be considerably smaller. Even for a system of uniformly banded 
equations such as occurs with the composite laminate fan blade shown in Figure 
4, a substantial reduction in computing time and some improvement in 
core-storage requirements are obtained with the profile solver package. 

For large problems such as the turbine blade model with platform as shown in 
Figure 2, the speed gain and reduction in memory needed make an in-core 
solution possible on a small computer. 

In a buckling analysis of the cylinder whose mesh is shown in Figure 3, the 
reduction in solution time is an order of magnitude more significant than the 
saving in storage requirements. 

Note that the tying option for the profile solver is implemented in a global 
manner. The algebraic constraint equations are embedded in the profile-form 
stiffness matrix after all the element matrices are assembled. On the other 
hand, the coordinate transformation option for the nodal degrees -of-freedom is 
handled on an element-by-element basis. Immediately after an element matrix is 
assembled, the entries subjected to coordinate transformations are manipulated 
before substituting into the global array. 
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Table 4 Performance Tests for the Profile Solver Solution Package 
in the MHOST Code 


Problem 1 SSME HPFTP Blade Model with 1025 Solid Elements (Type 7) and 
1575 Nodes (Mesh Shown in Figure 2) 

Profile Solver Band Solver 

Working Space Required 3,483,811 4,459,647 

Solution Time (seconds) 1466.594 NA* 

*Too large to be able to run on PRIME 9955 at MARC. 


Problem 2 Buckling Analysis of a Cylinder with 160 Shell Elements (Type 
75) and 176 Nodes (Mesh Shown in Figure 3) 

Profile Solver Band Solver 

Working Space 

- Static 

- Eigenvalue 

Solution Time (seconds) 

- Static 62.057 171.430 

- Eigenvalue 66.788 187.551 


Problem 3 Modal Analysis of Composite Laminate Fan Blade with 240 Shell 
Elements (Type 75) and 279 Nodes (Mesh Shown in Figure 4) 

Profile Solver Band Solver 

Working Space 824,341 1,054,259 

Solution Time (seconds) 24.794 93.764 


498,873 596,703 

1,044,773 1,340,375 
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Figure 2 




TEST CASE USING MHOST - PROFILE SOLVER; SSME HPFTP BLADE 


Finite Element Mesh 
PI atform) 


for Problem 1 (Turbine Blade Model with 



Figure 3 


Finite Element Mesh for Problem 2 (Buckling Analysis of a 
Cylinder) 



Figure 4 Finite Element Mesh for Problem 3 (Composite Laminate Fan Blade) 


4.4.2 Vibration Analysis of a Rotating Beam 

The stress stiffening and centrifugal mass options were validated by analyzing 
a cantilever beam rotating at a given angular velocity, with the geometry and 
loading conditions illustrated in Figure 5. Two finite element models were 
constructed on MHOST, one using the shell element (type 75) and one using the 
three-dimensional solid element (type 7). A reference solution was also 
established: a solution on NASTRAN using its shell element (QUAD4) was 
compared with a solution produced by MARC via its shell element (type 75). 

Both commercial codes produce almost identical vibration frequencies which 
establish the reference solution. As the values in Table 5 indicate, both 
MHOST solutions are almost identical to the reference solution. This also 
indicates that the selectivity integrated linear solid element is capable of 
accurately representing bending mode response. 
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Figure 5 Validation Model for Centrifugal Stiffening 

Table 5 Vibration Frequency-Cantilever Beam Under Centrifugal Loading 
(Column A: Without Stiffening; Column B: With Stiffening and 
Centrifugal Mass) 

NASTRAN MHOST 


Mode 

Shell 

El ement** 

Shell 

Element 

Solid El 

ement 


A 

B 

A 

B 

A 

B 

1* 

145.4 

227.6 

145.4 

228.7 

145.8 

229.9 

2 

618.9 

662.3 

628.6 

668.2 

595.8 

669.3 

3 

394.8 

1042.8 

924.5 

1077.0 

929.9 

1089.0 

4 

1115.0 

1199.0 

1316.7 

1 375.8 

1205.1 

1229.2 


* MARC shell element type 75 gives 146.2 Hz without centrifugal 
stiffening and 228.3 Hz with the stiffening and centrifugal mass 
effects included. 


** Element type QUAD4. 
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4.4.3 Large Displacement Analysis 

A plane stress test problem involving stretching an elastic-plastic plate is 
analyzed using the large displacement analysis capabilities of the MHOST code. 
Listing 1 shows the input data used for the calculations. The original finite 
element mesh, as well as a deformed configuration, are plotted in Figure 6. 

The solution obtained is smooth without any indication of numerical 
instabilities. Equivalent plastic strain and Mises equivalent stress plots are 
shown in Figures 7 and 8 respectively. 


Listing 1 Input Data for a Large Displacement Analysis 
by Plane Stress Elements 


ELASTIC PLASTIC ALUMINUM PLATE-8X8 PLANE STRESS-L 313 
*P0ST 
*N0DES 81 
♦ELEMENTS 64 
3 


♦CONSTITUTIVE 2 
♦HARD 2 
♦BOUNDARY 40 
*LARGE_DISPLACEMENT 
*FINITE_STRAIN 
♦STIFFENING 
♦LOUBIGNAC 3 1 3 
♦END 

♦WORKHARD 2 


45000 

.0 0.0000 

0.0000 


850000 

.0 10.000 

0.0000 


♦ITERATIONS 



20 

0.050 



♦COORDINATES 



1 

0.0000 

0.0000 

0.0000 

2 

0.0000 

1.0000 

0. 0000 

3 

0.0000 

2.0000 

0.0000 

4 

0.0000 

3.0000 

0. 0000 

5 

0.0000 

4.0000 

0.0000 


0.05000 
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6 

0.0000 

5.0000 

7 

0.0000 

6.0000 

8 

0.0000 

7.0000 

9 

0.0000 

8.0000 

10 

1.0000 

0.0000 

11 

1.0000 

1.0000 

12 

1.0000 

2.0000 

13 

1.0000 

3.0000 

14 

1.0000 

4.0000 

15 

1.0000 

5.0000 

16 

1.0000 

6.0000 

17 

1.0000 

7.0000 

18 

1.0000 

8.0000 

19 

2.0000 

0.0000 

20 

2.0000 

1.0000 

21 

2.0000 

2.0000 

22 

2.0000 

3.0000 

23 

2.0000 

4.0000 

24 

2.0000 

5.0000 

25 

2.0000 

6.0000 

26 

2.0000 

7.0000 

27 

2.0000 

8.0000 

28 

3.0000 

0.0000 

29 

3.0000 

1.0000 

30 

3.0000 

2.0000 

31 

3.0000 

3.0000 

32 

3.0000 

4.0000 

33 

3.0000 

5.0000 

34 

3.0000 

6.0000 

35 

3.0000 

7.0000 

36 

3.0000 

8.0000 

37 

4.0000 

0.0000 

38 

4.0000 

1.0000 

39 

4.0000 

2.0000 

40 

4.0000 

3.0000 


0.0000 
0.0000 
0.0000 
0.0000 
0. 0000 
0. 0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0. 0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
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41 

4.0000 

4.0000 

0.0000 

42 

4.0000 

5.0000 

0.0000 

43 

4.0000 

6.0000 

0.0000 

44 

4.0000 

7.0000 

0.0000 

45 

4.0000 

8.0000 

0.0000 

46 

5.0000 

0.0000 

0.0000 

47 

5.0000 

1.0000 

0.0000 

48 

5.0000 

2.0000 

0.0000 

49 

5.0000 

3.0000 

0.0000 

50 

5.0000 

4.0000 

0.0000 

51 

5.0000 

5.0000 

0.0000 

52 

5.0000 

6.0000 

0. 0000 

53 

5.0000 

7.0000 

0.0000 

54 

5.0000 

8.0000 

0.0000 

55 

6.0000 

0.0000 

0. 0000 

56 

6.0000 

1.0000 

0.0000 

57 

6.0000 

2.0000 

0. 0000 

58 

6.0000 

3.0000 

0.0000 

59 

6.0000 

4.0000 

0.0000 

60 

6.0000 

5.0000 

0.0000 

61 

6.0000 

6.0000 

0.0000 

62 

6.0000 

7.0000 

0.0000 

63 

6.0000 

8.0000 

0.0000 

64 

7.0000 

0. 0000 

0.0000 

65 

7.0000 

1.0000 

0.0000 

66 

7.0000 

2.0000 

0.0000 

67 

7.0000 

3.0000 

0.0000 

68 

7.0000 

4.0000 

0.0000 

69 

7.0000 

5.0000 

0.0000 

70 

7.0000 

6.0000 

0.0000 
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71 

7.0000 

72 

7.0000 

73 

8.0000 

74 

8.0000 

75 

8.0000 

76 

8.0000 

77 

8.0000 

78 

8.0000 

79 

8.0000 

80 

8.0000 

81 

8.0000 

♦ELEMENTS 3 

1 

1 10 

2 

2 11 

3 

3 12 

4 

4 13 

5 

5 14 

6 

6 15 

7 

7 16 

8 

8 17 

9 

10 19 

10 

11 20 

11 

12 21 

12 

13 22 

13 

14 23 

14 

15 24 

15 

16 25 

16 

17 26 

17 

19 28 

18 

20 29 

19 

21 30 

20 

22 31 


7.0000 0.0000 

8.0000 0.0000 

0.0000 0.0000 

1.0000 0.0000 

2.0000 0.0000 

3.0000 0.0000 

4.0000 0.0000 

5.0000 0.0000 

6.0000 0.0000 

7.0000 0.0000 

8.0000 0.0000 


11 

2 

12 

3 

13 

4 

14 

5 

15 

6 

16 

7 

17 

8 

18 

9 

20 

11 

21 

12 

22 

13 

23 

14 

24 

15 

25 

16 

26 

17 

27 

18 

29 

20 

30 

21 

31 

22 

32 

23 


45 


21 

23 

32 

33 

24 

22 

24 

33 

34 

25 

23 

25 

34 

35 

26 

24 

26 

35 

36 

27 

25 

28 

37 

38 

29 

26 

29 

38 

39 

30 

27 

30 

39 

40 

31 

28 

31 

40 

41 

32 

29 

32 

41 

42 

33 

30 

33 

42 

43 

34 

31 

34 

43 

44 

35 

32 

35 

44 

45 

36 

33 

37 

46 

47 

38 

34 

38 

47 

48 

39 

35 

39 

48 

49 

40 

36 

40 

49 

50 

41 

37 

41 

50 

51 

42 

38 

42 

51 

52 

43 

39 

43 

52 

53 

44 

40 

44 

53 

54 

45 

41 

46 

55 

56 

47 

42 

47 

56 

57 

48 

43 

48 

57 

58 

49 

44 

49 

58 

59 

50 

45 

50 

59 

60 

51 

46 

51 

60 

61 

52 

47 

52 

61 

62 

53 

48 

53 

62 

63 

54 

49 

55 

64 

65 

56 

50 

56 

65 

66 

57 


! 

! 


46 



51 

57 

66 

67 

58 

52 

58 

67 

68 

59 

53 

59 

68 

69 

60 

54 

60 

69 

70 

61 

55 

61 

70 

71 

62 

56 

62 

71 

72 

63 

57 

64 

73 

74 

65 

58 

65 

74 

75 

66 

59 

66 

75 

76 

67 

60 

67 

76 

77 

68 

61 

68 

77 

78 

69 

62 

69 

78 

79 

70 

63 

70 

79 

80 

71 

64 

71 

80 

81 

72 


♦BOUNDARY 
1 1 
1 2 
2 1 

3 1 

4 1 

5 1 

6 1 

7 1 

8 1 

9 1 

10 2 
19 2 

28 2 
37 2 

46 2 
55 2 
64 2 

73 1 0.8000 

73 2 
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74 

1 

0.8000 

74 

2 


75 

1 

0.8000 

75 

2 


76 

1 

0.8000 

76 

2 


77 

1 

0.8000 

77 

2 


78 

1 

0.8000 

78 

2 


79 

1 

0.8000 

79 

2 


80 

1 

0.8000 

80 

2 


81 

1 

0.8000 

81 

2 



♦PROPERTIES 3 

C STANDARD ALUMINUM MATERIAL 
1 81 0.05 10.0E6 0.330 

♦PRINTOPTION 

INCREMENTALDISPLACEMENT NODE 


TOTALDISPLACEMENT NODE 

REACTION NODE 

PLASTIC NODE 

STRESS NODE 

STRAIN NODE 

FORCE NODE 

♦END 

♦AUTOINCREMENT 

9 

♦END 

♦STOP 
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FIGURE 7 ELASTIC-PLASTIC MEMBRANE: CONTOURS OF EQUIVALENT PLASTIC STRAIN 
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FIGURE 8 ELASTIC-PLASTIC MEMBRANE: CONTOURS OF MISES STRESS 


A disc upsetting problem involving both large displacement and large strain 
behavior is used as a second example. This is a standard problem solved 
repeatedly by MARC, and the reference solution is readily available. As the 
initial mesh and deformed configuration plots (Figure 9) indicate, a large 
amplitude hourglassing is observed. This is due to the nearly incompressible 
response of the material and the same phenomenon with less significant 
amplitude is observable in a standard displacement solution. Further 
investigation on stabilization of the displacement field is necessary to 
filter out this spurious mode. However, the stress field obtained with the 
present capability is stable and quantitatively agrees well with the solution 
obtained by MARC via the conventional displacement method (Figure 10). 

For comparison, the same disc upsetting problem is studied using the 
three-dimensional solid element (Figure 11). As indicated in the stress 
contour plot (Figure 11) at an early stage of the deformation process, the 
model is capable of reproducing an axi symmetric stress field quantitatively 
similar to that generated using the axi symmetric element. 

As an example of MHOST's large displacement analysis capability with shell 
elements, a one-quarter model of a clamped square plate under uniform pressure 
loading is studied. The results are plotted in Figure 12, together with 
solutions of the same problem available in the literature [WAY (1938); KAWAI, 
YOSHIMURA (1961); HUGHES, LIU (1981); and NAGTEGAAL, SLATER (1981)]. The 
comparison indicates that the mixed iterative implementation of shell elements 
results in more flexible behavior compared with the compatible model 
(Kawai-Yoshimura solution) and with models based on the Reissner-Mindl in 
theory (Hughes-Liu solution by uniform reduced integration and heterosis 
elements). Qualitatively, the behavior of the MHOST shell element is similar 
to the Nagtegaal -Slater solution via the SLICK element derived from the 
discrete Kirchhoff assumption (MARC element type 72). A visible difference is 
shown between the solution obtained with 'full' Newton iteration (in which the 
displacement preconditioner is repeatedly updated) and that obtained via the 
secant-Newton update with line search. This difference may be due to 
ill -conditioning of the preconditioner caused by transverse shear constraint 
terms. 
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Figure 9 



Compression of Solid Axi symmetric Cylinder; Deformed Shape 
(Hourglassing) 
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FIGURE 10 AXISYMMETRIC UPSETTING: EQUIVALENT MISES TENSILE STRESS-MHOST VS 
MARC RESULTS 







E = 2 x TO 4 kg/mm 2 
v = 0.3 
L = 1000 mm 
h = 2 mm 
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Figure 12 Load Deflection Curves of a Clamped Square Plate Under Uniform 
Pressure Loading 
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4.4.4 Application of the Subelement Method for Inelastic Analysis 

A plane stress problem of a square plate with a circular hole is used to 
demonstrate the inelastic solution capability of the MHOST subelement 
procedures. The global finite element mesh and the subelement mesh 
representing a hole embedded in the element at the center of a plate is 
illustrated in Figure 13. A global finite element model explicitly including 
the hole in the same way as the subelement grid is shown in Figure 14. 



SUBELEMENT MESH FOR AN EMBEDDED HOLE 


ELEMENT DIVIDED INTO A SUBELEMENT 



Figure 13 Global and Subelement Meshes for an Embedded Hole 
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TEST PROBLEM FOR A PLATE WITH A HOLE IN PLANE STRESS 


Figure 14 Global Finite Element Mesh 


The input data files for both the global element mesh and the subelement mesh 
models are given in Listings 2 and 3. Note that the amount of data to 
represent a hole is significantly less for the subelement mesh; this 
represents a major improvement in user friendliness. 
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Listing 2 Global Finite Element Model 

TEST PROBLEM FOR A PLATE WITH A HOLE ( PLANE STRESS ) 

*POST 

*NOECHO 


♦CONSTITUTIVE 0 


♦ELEMENTS 

28 




3 




♦NODES 


40 



*BOUNDARY_C 

19 



♦HARDENING 

5 



♦FORCE 


10 



♦END 





♦ITERATION 




20 


.01 

10 . 

0 

♦ELEMENT 

3 



1 

1 

5 

6 

2 

2 

2 

6 

21 

17 

3 

17 

21 

7 

3 

4 

3 

7 

8 

4 

5 

5 

18 

22 

6 

6 

6 

22 

27 

25 

7 

6 

25 

39 

21 

8 

21 

39 

37 

7 

9 

37 

35 

24 

7 

10 

7 

24 

20 

8 

11 

25 

27 

28 

26 

12 

25 

26 

40 

39 

13 

39 

40 

38 

37 

14 

38 

36 

35 

37 

15 

18 

9 

10 

22 

16 

22 

10 

29 

27 

17 

29 

10 

23 

31 

18 

31 

23 

11 

33 

19 

33 

11 

24 

35 
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20 

24 

11 

12 

20 

21 

27 

29 

30 

28 

22 

30 

29 

31 

32 

23 

32 

31 

33 

34 

24 

34 

33 

35 

36 

25 

9 

13 

14 

10 

26 

10 

14 

19 

23 

27 

23 

19 

15 

11 

28 

11 

15 

16 

12 


♦COORDINATE 

1 0.0 0.0 

2 0.0 1.0 

3 0.0 2.0 

4 0.0 3.0 

5 1.0 0.0 

6 1.0 1.0 

7 1.0 2.0 

8 1.0 3.0 

9 2.0 0.0 

10 2.0 1.0 

11 2.0 2.0 

12 2.0 3.0 

13 3.0 0.0 

14 3.0 1.0 

15 3.0 2.0 

16 3.0 3.0 

17 0.0 1.5 

18 1.5 0.0 

19 3.0 1.5 

20 1.5 3.0 

21 1.0 1.5 

22 1.5 1.0 

23 2.0 1.5 

24 1.5 2.0 
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25 

1.2146 

1.21 46 


26 

1.4293 

1 .4293 


27 

1.5000 

1 .2000 


28 

1 . 5000 

1.4000 


29 

1.7854 

1.21 46 


30 

1.5707 

1.4293 


31 

1.8000 

1 . 5000 


32 

1.6000 

1.5000 


33 

1.7854 

1 . 7854 


34 

1.5707 

1.5707 


35 

1.5000 

1.8000 


36 

1 . 5000 

1.6000 


37 

1.21 46 

1 . 7854 


38 

1.4293 

1.5707 


39 

1 .2000 

1 . 5000 


40 

1.4000 

1.5000 

♦PROPERTY 



1 ■ 

40 

1.0 1 

*BOUNDARY_C 
1 1 

1 2 

2 1 

3 1 

4 1 

17 1 


C 

13 

1 

0.01 

C 

14 

1 

0.01 

C 

15 

1 

0.01 

C 

16 

1 

0.01 


. 0E6 0.3 
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*FORCE 


13 

1 

10000.0 

14 

1 

15000.0 

19 

1 

1 0000. 0 

15 

1 

15000.0 

16 

1 

10000.0 


*PRINT 

TOTAL 

STRESS 

REACTION 

*END 

*STOP 


Listing 3 Subelement Model 

TEST PROBLEM FOR AN EMBEDDED HOLE ( PLANE STRESS ) 
♦SECA 

*LOUB 3 1 3 

♦NOECHO 

*EMBED 

♦CONSTITUTIVE 0 

♦ELEMENTS 9 

3 

C 101 

♦NODES 1 6 

*BOUNDARY_C 9 

♦FORCE 4 

♦END 

♦ITERATION 

20 .1 10.0 1.00 
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♦ELEMENT 
1 1 

2 2 

3 3 

4 5 

5 6 

6 7 

7 9 

8 10 

9 11 

♦COORDINATE 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

♦PROPERTY 
1 16 

*BOUNDARY_C 
1 1 

1 2 

2 1 

3 1 

4 1 


3 


5 

6 

2 

6 

7 

3 

7 

8 

4 

9 

10 

6 

10 

11 

7 

11 

12 

8 

13 

14 

10 

14 

15 

11 

15 

16 

12 

0.0 


0.0 

O 

• 

o 


1.0 

0.0 


2.0 

0.0 


3.0 

1.0 


O 

• 

o 

1.0 


1.0 

1.0 


2.0 

1.0 


3.0 

2.0 


0.0 

2.0 


1.0 

2.0 


2.0 

2.0 


3.0 

3.0 


0.0 

3.0 


1.0 

3.0 


2.0 

3.0 


3.0 


1.0 


. 0E6 0.3 


68 



c 

13 

1 

0.01 

c 

14 

1 

0.01 

c 

15 

1 

0.01 

c 

16 

1 

0.01 


♦FORCE 

13 1 10000.0 

14 1 20000.0 

15 1 20000.0 

16 1 10000.0 

♦HOLE 

5 3 0.2 

♦PRINT 

TOTAL 

STRESS 

REACTION 

♦END 

♦STOP 


The elastic stress field obtained by the global solution is shown in Figure 15 
and qualitatively agrees with the classical solution found in many textbooks, 
e.g., SAVIN (1951). The elastic solution after one iteration, as summarized in 
Table 6, indicates that the quadratic subelement produces the most accurate 
results for the stress concentration factor. 


69 


6.092 + 3 


ORIGINAL PAGE 
COLOR PHOTOGRAPH 



CO 











+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

00 



CO 

CD 

in 

T— 

CO 


o 

CO 

LD 

o 


CD 

CO 

CO 

CO 


CM 

r^ 

*— 

in 

CO 

CD 

O) 

CO 

CO 

o 

CO 

r>- 

o 


CT> 


r* 

r- 

c\i 

csi 

CO 

CO 

CO 




□ 


PRECEDING PAGE BLANK NOT FILMED 




71 


FIGURE 1 5 ELASTIC STRESS DISTRIBUTION FOR A PLATE WITH A HOLE 



Table 6 Nodal Stress Concentration Factor at the Edge 
of a Circular Hole After One Iteration 



Global Solution 

Subelement Solution 

Point A 


Linear 

Quadratic 

°x 

2.2615 

2.31 98 

2.7087 

°y 

0.5916 

0.5977 

0.3929 

T xy 

0.0002 

0.0000 

0.0003 

Point B 





0.1752 

0.4339 

-0.2702 

°y 

-0.0813 

-0.2995 

-0.5958 

T xy 

0.0000 

0.0000 

0.0000 


For the elastic-plastic analyses, the material is assumed to be bilinear and 
is defined by *W0RKL0AD data as: 

Stress Plastic Strain 

2.0 0.0 

4.0 0.4 

Note that the stress here is normalized by the uniform horizontal stress for 
the solution of the square plate without a hole. As is well-known from 
standard elasticity, amplification of the horizontal stress value by a factor 
of 3.0 is expected at Point A under pure elastic behavior. 

Both linear and quadratic subelement results, as well as global finite element 
results, with the same mesh, are compared. The monotonic convergence 
characteristics of the subelement method can be seen in these comparisons, as 
summarized in Table 7. 
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Table 7 Elastic-Plastic Analysis of a Square Plate with a Circular Hole; 
Nodal Stress at the Edge of the Circular Hole at Maximum SCF 
Location (Point A) 

Linear Quadratic 


Global Solution 

Subel ement 

Solution 

Subelement Solution 

Elastic Elastic- 
Plastic 

Elastic 

Elastic- 

Plastic 

Elastic Elastic- 

Plastic 


Equivalent 
Plastic Strain 

0.0 

0.01 91 

0.0 

0.0320 

0.0 

0.0513 

Stress 

crxx 

<ryy 

crxy 

2.3492 
0.31 60 
0.0034 

2.2786 

0.4339 

0.0049 

2.5064 

0.4327 

0.0010 

2.41 51 
0.6682 
0.0009 

2.7087 

0.3929 

0.0000 

2.4624 

0.4929 

0.0000 


Note here that iterations were carried out until a convergence tolerance of 
10% in terms of relative residual was met. This solution procedure improved 
the elastic stress results compared to elastic results obtained earlier. The 
same convergence criteria were used for all calculations shown in Table 7. In 
comparison with solutions reported previously, the quadratic subelement method 
exhibited superior performance, at least for this particular example of an 
embedded hole. 


74 



5.0 REFERENCES 


Babuska, I., 0. C. Zienkiewicz, J. Gago and E. R. de A. Oliveira (eds.). 
Accuracy Estimates and Adaptive Refinements in Finite Element Computations , 
Wiley, Chichester, 1986. 

Babuska, I., "The Finite Element Method with Lagrange Multipliers," Numer. 
Math. , 20, 1973, pp. 179-192. 

Bathe, K. J. and E. L. Wilson, Numerical Methods in Finite Element Analysis , 
Prentice Hall, Englewood Cliffs, New Jersey, 1976. 

Bazeley, G. P., Y. K. Cheung, B. M. Irons and 0. C. Zienkiewicz, "Triangular 
Elements in Plate Bending. Conforming and Non-Conforming Solutions," Proc. 1st 
Conf. Matrix Methods in Structural Mechanics, AFFDL-TR-CC-80, 1966, pp. 
547-576. 

Belytschko, T. , "A Review of Recent Developments in Plate and Shell Elements," 
Computational Mechanics - Advances and Trends (ed. A. K. Noor), ASME AMD - 
V'ol. 75, V986, pp'.~21 7-232. 

Benson, D. J. and J. 0. Hallquist, "A Simple Rigid Body Algorithm for 
Structural Dynamics Programs," Int. J. Num. Meth. Eng. , 22, 1986, pp. 723-749. 

Brezzi , F., "On the Existence, Uniqueness and Approximation of Saddle-Point 
Problems Arising from Lagrange Multipliers," RAIRO , 22, 1974, pp. 129-151. 

Brown, P. N. and A. C. Hindmarsh, "Matrix-Free Methods for Stiff System of 
ODE's," SIAM J. Numer. Anal. , 23, 1986, pp. 610-638. 

Fyhre, D. and T. J. R. Hughes, "A Literature Survey on Finite Element 
Methodology for Turbine Engine Hot Section Components," MARC Analysis Research 
Corporation, Palo Alto, California, 1983. 

Gaudreau, G. L., D. J. Benson, J. 0. Hallquist, G. J. Kay, R. W. Rosinsky and 
S. J. Sachett, "Supercomputers for Engineering Analysis," Computational 
Mechanics - Advances and Trends (ed. A. K. Noor), ASME AMD - Vol . 75, 1986, 
pp. 117 - 131 . 

Hayes, L. J. and P. Devloo, "A Vectorized Version of a Spare Matrix-Vector 
Multiply," Int. J. Num. Meth. Eng. , 23, 1986, pp. 1043-1056. 

Hu, H. C., "On Some Variational Principles in the Theory of Elasticity and 
Plasticity," Scintia Sinica , 4, 1955, pp. 33-54. 

Hughes, T. J. R. and W. K. Liu, "Nonlinear Finite Element Analysis of Shells; 
Part I: Three Dimensional Shells," Comp. Meth. Appl . Mech. Eng., 26, 1981, pp. 
331 -362. 


75 


Hughes, T. J. R. , "Analysis of Transient Algorithms with Particular Reference 
to Stability Behavior," Computational Methods for Transient Analysis (eds. T. 
Belytschko and T. J. R. Hughes), Chap. 2, North-Holland, Amsterdam, 1983. 

Hughes, T. J. R. and F. Shakie, "Pseudo-Corner Theory: A Simple Enhancement of 
J 2 - Flow Theory for Applications Involving Non-Proportional Loading," Eng. 
Comp ♦ , 3, 1986, pp. 116-120. 

Hughes, T. J. R., R. M. Ferencz and J. 0. Hallquist, "Large-Scale Vectorized 
Implicit Calculations in Solid Mechanics on a CRAY X-MP/48 Utilizing EBE 
Preconditioned Conjugate Gradients," Computational Mechanics - Advances and 
Trends (ed. A. K. Noor), ASME AMD - Vol. 75, 1986, pp. 233-279. 

Jameson, A., "Current Status and Future Directions of Computational 
Transonics," Computational Mechanics - Advances and Trends (ed. A. K. Noor), 
ASME AMD - Vol. 76, 1986, pp. 329-368. 

Kawai , T. and N. Yoshimura, "Analysis of Large Deflection of Plates by Finite 
Element Method," Int, J. Num. Meth. Eng. , 1 , 1969, pp. 123-133. 

Le Tallec, P., "Mixed Numerical Methods in Viscoplasticity," Proc. 1st WCCM , 
Vol. 1, Texas Institute for Computational Mechanics, 1986. 

Lohner, R. , "FEM-FCT and Adaptive Refinement Schemes for Strongly Unsteady 
Flows," Numerical Methods for Compressible Flows - Finite Difference, Element 
and Volume Techniques (eds. T. E. Tezduyar and T. J. R. Hughes), ASME AMD - 
Vol .Tff.H 986, pp. 93-114. 

Nagtegaal , J. C. and J. G. Slater, "A Simple Non-Compatible Thin Shell Element 
Based on Discrete Kirchhoff Theory," Non! inear Finite Element Analysis of 
Plates and Shells (ed. T. J. R. Hughes, A. Pifko and A. Jay), ASME AMD - Vol. 
48, 1981 , pp. 167-1 92 . 

Nagtegaal, J. C., S. Nakazawa and M. Tateishi, "On the Construction of Optimal 
Mindlin Type Shell Elements," Finite Element Methods for Plate and Shell 
Structures, Vol. 1: Element TecHnology (eds. T. J. R. Hughes and E. Hinton), 
Pineridge Press, Swansea, 1986, pp. 348-364. 

Nakazawa, S., 3-D Inelastic Analysis Methods for Hot Section Components, Third 
Annual Status Report; Volume 1: Special Finite Element Models, NASA CR-179494, 
Pratt & Whitney PWA-5940-46, 1986. 

Nakazawa, S. and J. C. Nagtegaal, "An Efficient Numerical Procedure for 
Nonlinear Analysis of Turbine Engine Hot Section Components," presented at 3rd 
Symp. on Nonlinear Constitutive Relations for High Temperature Applications , 
Akron, Ohio, 1986. 

Nakazawa, S. and M. S. Spiegel, "Mixed Finite Elements and Iterative Solution 
Processes," presented at 1st WCCM , Austin, Texas, 1986. 

Needleman, A., "Finite Element Analysis of Failure Modes in Ductile Solids," 
Proc. 1st WCCM, Vol. 1, Texas Institute for Computational Mechanics, 1986. 


76 


Nour-Omid, B. and M. Ortiz, "Concurrent Solution Methods for Finite Element 
Equations," Proc. 1st WCCM , Vol . 1, Texas Institute for Computational 
Mechanics, 1 986. 

Ortiz, M. and Y. Leroy, "A Finite Element Method for Localized Failure 
Analysis," Proc, 1st WCCM , Vol. 1, Texas Institute for Computational 
Mechanics, 1 986. 

Ortiz, M. and J. C. Simo, "An Analysis of a New Class of Integration 
Algorithms for Elastoplastic Constitutive Relations," Int. J. Num. Meth. Eng., 
23, 1986, pp. 355-366. 

Peraire, J., K. Morgan and 0. C. Zienkiewicz, "Convection Dominated Problems," 
Numerical Methods for Compressible Flows - Finite Difference, Element and 
Volume Techniques (eds. T. E. Tezduyar and T. J. R. Hughes), ASME AMD - Vol . 
W, 1986, pp. 1 29-147. 

Pian, T. H. H. (ed.). Special Memorial Issue Dedicated to Prof. Kyuichiro 
Washizu, Computers and Structures , 19 (1/2), 1984. 

Rawtani , S. and M. A. Dokainish, "Natural Frequencies of Rotating Low Aspect 
Ratio Turbomachinery Blades," AIAA Paper 71-374, Anaheim, California, 1971. 

Savin, G. N., Stress Concentration Around Holes (in Russian, 1951), English 
Edition by Pergamon Press, London, 1961. 

Simo, J. C., "Regular Perturbation Schemes for the Iterative Solution of 
Symmetric Eigenvalue Problems," Manuscript, Applied Mechanics Division, 
Stanford University, 1985. 

Simo, J. C., "A Framework for Finite Strain Elastoplasticity Based on Maximum 
Plastic Dissipation and the Multiplicative Decomposition. Part I: Continuum 
Formulation; Part II: Computational Aspects," Manuscript, Applied Mechanics 
Division, Stanford University, 1985. 

Stummel , F. , "The Limitations of Patch Test," Int. J. Num. Meth. Eng., 15, 
1980, pp. 177-188. 

Taylor, R. L., "Solution of Linear Equations by a Profile Solver," Eng. Comp., 
2, 1985, pp. 344-350. 

Taylor, R. L., J. C. Simo, 0. C. Zienkiewicz and A. C. H. Chen, "The Patch 
Test - A Condition for Assessing FEM Convergence," Int. J. Num. Meth. Eng., 

22, 1986, pp. 39-62. 

Thomas, J. and C. A. Mota Soares, "Finite Element Analysis of Rotating 
Shells," ASME Paper 73-DET-94, 1973. 

Utku, S., M. Salama and R. J. Melosh, "Concurrent Cholesky Factorization of 
Positive Definite Banded Hermitian Matrices," Int. J. Num. Meth. Eng., 23, 
1986, pp. 21 37-2152. 


77 


Utku, S., R. J. Melosh, M. Salama and H. T. Chang, "Problem Decomposition for 
Parallel Processing," Proc. 1st WCCM , Vol. 1, Texas Institute for 
Computational Mechanics, 1986. 

Washizu, K., On the Variational Principles of Elasticity and Plasticity , 
Technical Report 25-18, Aeroelastic and Struct. Res. Lab., Massachusetts 
Institute of Technology, 1955. 

Washizu, K., Variational Methods in Elasticity and Plasticity , Pergamon Press, 
Oxford, 1974. 

Way, S., "Uniformly Loaded, Clamped, Rectangular Plates with Large 
Deformation," Proc. 5th Int'l. Cong. Appl . Mech., Cambridge, Massachusetts, 
1938. 

William, K., N. Sobh and S. Sture, "Bifurcation Studies of Elastic-Plastic 
Solids," Proc. 1st WCCM, Vol. 1, Texas Institute for Computational Mechanics, 

1 986. 

Wilson, R. B., M. J. Bak, S. Nakazawa and P. K. Banerjee, 3-D Inelastic 
Analysis Methods for Hot Section Components (Base Program), First Annual 
Status "Report' , NA5ATIT-1T47(flJr 1 984: 

Wilson, R. B., M. J. Bak, S. Nakazawa and P. K. Banerjee, 3-D Inelastic 
Analysis Methods for Hot Section Components (Base Program), Second Annual 
Status Report , NASA Cft-1 75060, 1985. 

Zienkiewicz, 0. C., The Finite Element Method, 3rd Edition, McGraw-Hill, New 
York, 1977. 

Zienkiewicz, 0. C. and S. Nakazawa, "On Variational Formulation and Its 
Modifications for Numerical Solution," Computers and Structures, 19, 1984, pp. 
303-313. 

Zienkiewicz, 0. C., X. K. Li and S. Nakazawa, "Dynamic Transient Analysis by a 
Mixed, Iterative Method," Int. J. Num. Meth. Eng. , 23, 1986, pp. 1343-1353. 

Zienkiewicz, 0. C., Three-Field Mixed Approximation and the Plate Bending 
Problem, Report No. C/R/566/86,' Institute for Numerical Method in Engineering, 
University College of Swansea, 1986. 

Zienkiewicz, 0. C., S. Qu, R. L. Taylor and S. Nakazawa, "The Patch Test for 
Mixed Formulations," Int. J. Num. Meth. Eng. , 23, 1986, pp. 1873-1883. 

Zienkiewicz, 0. C. and J. Z. Zhu, "A Simple Error Estimator and Adaptive 
Procedure for Practical Engineering Analysis," Int. J. Num. Meth. Eng. , 24, 

1987, pp. 337-357. 

Zyvoloski, G. , "Incomplete Factorization for Finite Element Methods," Int. J. 
Num. Meth. Eng., 23, 1986, pp. 1101-1109. 


78 



DISTRIBUTION LIST 


NASA-Lewis Research Center 
Attn: Contracting Officer, MS 500-313 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: Tech Rpt Control Office, MS 60-1 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: Tech Utilization Office, MS 7-3 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: AFSC Liaison Office, MS 501-3 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: Div. Contract File, MS 49-6 
21000 Brookpark Road 
Cleveland, OH 44135 (2 copies) 

NASA-Lewis Research Center 
Attn: Library, MS 60-3 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: C. C. Chamis, MS 49-8 
21 000 Brookpark Road 
Cleveland, OH 44135 (2 copies) 

NASA-Lewis Research Center 
Attn: L. Berke, MS 49-6 
21 000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: R. E. Kielb, MS 49-8 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: G. R. Halford, MS 49-7 
21000 Brookpark Road 
Cleveland, OH 44135 


NASA-Lewis Research Center 
Attn; D. A. Hopkins, MS 49-8 
21000 Brookpark Road 
Cleveland, OH 44135 (2 copies) 

NASA-Lewis Research Center 
Attn: R. H. Johns, MS 49-8 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: D. E. Sokolowski, MS 49-7 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: R. L. Thompson, MS 49-8 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA-Lewis Research Center 
Attn: R. C. Bill, MS 77-12 
21000 Brookpark Road 
Cleveland, OH 44135 

NASA Headquarters 
Attn: NHS-22/Library 
600 Independence Avenue, S.W. 
Washington, DC 20546 

NASA Headquarters 
Attn: R/S. L. Venneri 
600 Independence Avenue, S.W. 
Washington, DC 20546 

NASA Headquarters 
Attn: R/M. S. Hirschbein 
600 Independence Avenue, S.W. 
Washington, DC 20546 

NASA S&T Information Facility 
Attn: Acquisition Dept. (10 copies) 
P. 0. Box 8757 

Balt. -Wash. Int. Airport, MD 21240 

NASA Ames Research Center 
Attn: Library, MS 202-3 
Moffett Field, CA 94035 


79 


DISTRIBUTION LIST (continued) 


NASA Goddard Space Flight Center 
Attn: 252/Library 
Greenbelt, MD 20771 

NASA John F. Kennedy Space Center 
Attn: Library, AD-CSO-1 
Kennedy Space Center, FL 32931 


NASA Langley Research Center 
Attn: Library, MS 185 
Hampton, VA 23665 


NASA Langley Research Center 
Attn: J. Starnes 
Hampton , VA 23665 


NASA Lyndon B. Johnson Space Center 
Attn: JM6/Library 
Houston, TX 77001 


NASA George C. Marshall Space 
Flight Center 
Attn: L. Kief ling 

Marshall Space Fit. Center, AL 35812 

Jet Propulsion Laboratory 
Attn: B. Wada 
4800 Oak Grove Drive 
Pasadena, CA #91103 

Air Force Aeronautical Propulsion 
Laboratory 
Attn: Z. Gershon 
Wright-Patterson AFB, OH 45433 

Air Force Systems Command 
Attn: C. W. Cowie 
Aeronautical Systems Division 
Wright-Patterson AFB, OH 45433 

Air Force Systems Command 
Attn: R. J. Hill 
Aeronautical Systems Division 
Wright-Patterson AFB, OH 45433 


U.S. Army Ballistics Research Lab. 
Attn: Dr. D. F. Haskell, DRXBR-BM 
Aberdeen Proving Ground, MD 21005 

Mechanics Research Laboratory 
Attn: Dr. E. M. Lenoe 
Army Mat'ls 4 Mech. Research Center 
Watertown, MA 02172 

National Bureau of Standards 
Attn: R. Mitchell 
Engineering Mechanics Section 
Washington, DC 20234 

U.S. Army Missile Command 
Attn: Document Section 
Redstone Scientific Info Center 
Redstone Arsenal , AL 35808 

Commanding Officer 
U.S. Army Research Office (Durham) 
Box CM, Duke Station (Attn: Library) 
Durham, NC 27706 

Department of the Army 
Attn: AMCRD-RC 
U.S. Army Material Command 
Washington, DC 20315 

Bureau of Naval Weapons 
Department of the Navy 
Attn: RRRE-6 
Washington, DC 20360 

Commander, U.S. Naval Ordnance Lab. 
Attn: Library 
White Oak 

Silver Springs, MD 20910 

Director, Code 6180 
U.S. Naval Research Laboratory 
Attn: Library 
Washington, DC 20390 

Naval Air Propulsion Test Center 
Aeronautical Engine Department 
Attn: Mr. James Sal vino 
Trenton, NJ 08628 


80 



DISTRIBUTION LIST (continued) 


Denver Federal Center 
U.S. Bureau of Reclamation 
Attn: P. M. Lorenz 
P. 0. Box 25007 
Denver, CO 80225 

FAA, ARD-520 

Attn: Commander John J. Shea 
2100 Second Street, SW 
Washington, DC 20591 

Federal Aviation Administration DOT 
Office of Aviation Safety, FOB 1 OA 
Attn: Mr. John H. Enders 
800 Independence Avenue, SW 
Washington, DC 20594 

Federal Aviation Administration 
Code ANE-214, Propulsion Section 
Attn: Mr. Robert Berman 
12 New England Executive Park 
Burlington, MA 01803 

National Transportation Safety Board 
Attn: Edward P. Wizniak, MS TE-20 
800 Independence Avenue, SW 
Washington, DC 20594 

Aeronautical Research Association 
of Princeton, Inc. 

Attn: Dr. Thomas McDonough 
P. 0. Box 2229 
Princeton, NJ 08540 

Aerospace Corporation 
Attn: Library-Documents 
2400 E. El Segundo Blvd. 

Los Angeles, CA 90045 

Allison Gas Turbine Operation 
General Motors Corporation 
Attn: Mr. William Springer 
Speed Code T3, Box 894 
Indianapolis, IN 46206 


Allison Gas Turbine Operation 
General Motors Corporation 
Attn: Mr. L. Snyder 
Speed Code T3, Box 894 
Indianapolis, IN 46206 

AVCO Lycoming Division 
Attn: Mr. Herbert Kaehler 
550 South Main Street 
Stratford, CT 06497 

Boeing Aerospace Company 
Impact Mechanics Lab 
Attn: Dr. R. J. Bristow 
P. 0. Box 3999 
Seattle, WA 98124 

Boeing Commercial Airplane Company 
Attn: Dr. Ralph B. McCormick 
P. 0. Box 3707 
Seattle, WA 98124 


Bell Aerospace 
Attn: R. A. Gellatly 
P. 0. Box 1 
Buffalo, NY 14240 

Bell Aerospace 
Attn: S. Gellin 
P. 0. Box 1 
Buffalo, NY 14240 


Beech Aircraft Corporation, 

Plant 1 

Attn: Mr. M. K. O'Connor 
Wichita, KA 67201 

Boeing Commercial Airplane Company 
Attn: David T. Powell, MS 73-01 
P. 0. Box 3707 
Seattle, WA 98124 


81 


DISTRIBUTION LIST (continued) 


Boeing Commercial Airplane Company 
Attn: Dr. John H. Gerstle 
P. 0. Box 3707 
Seattle, WA 98124 

Boeing Company 

Attn: Mr. C. F. Tiffany 

Wichita, KA 67201 

Douglas Aircraft Company 

Attn: M. A. O'Connor, Jr., MS 36-41 

3855 Lakewood Blvd. 

Long Beach, CA 90846 

Garrett Ai Research Manufacturing Co. 

Attn: L. A. Matsch 

111 S. 34th Street 

P. 0. Box 5217 

Phoenix, AZ 85010 

General Dynami cs 
Attn: Library 
P. 0. Box 748 
Fort Worth, TX 76101 

General Dynami cs/Convair Aerospace 
Attn: Library 
P. 0. Box 1128 
San Diego, CA 92112 

General Electric Company 
Attn: Dr. L. Beitch, MS K221 
Interstate 75, Bldg. 500 
Cincinnati, OH 45215 

General Electric Company 
Attn: Dr. R. L. Mcknight, MS K221 
Interstate 75, Bldg. 500 
Cincinnati, OH 45215 


General Electric Company 
Attn: Dr. M. Roberts, MS K221 
Interstate 75, Bldg. 500 
Cincinnati, OH 45215 


General Electric Company 
Aircraft Engine Group 
Attn: Mr. Herbert Garten 
Lynn, MA 01902 

General Motors Corporation 
Attn: R. J. Tippet 
Warren, MI 48090 

Grumman Aircraft Engineering 
Corporation 
Attn: Library 

Bethpage, Long Island, NY 11714 

Grumman Aircraft Engineering 
Corporation 
Attn: H. A. Armen 
Bethpage, Long Island, NY 11714 


Lawrence Livermore Laboratory 
Attn: Library 
P. 0. Box 808, L-421 
Livermore, CA 94550 

Lawrence Livermore Laboratory 
Attn: M. L. Wilkins 
P. 0. Box 808, L-421 
Livermore, CA 94550 

Arthur D. Little 
Acorn Park 
Attn: P. D. Hilton 
Cambridge, MA 02140 

Lockheed California Company 
Attn: Mr. D. T. PI and 
P. 0. Box 551 

Dept. 73-31, Bldg. 90, Plant A-l 
Burbank, CA 91520 

Lockheed California Company 
Attn: Mr. Jack E. Wignot 
P. 0. Box 551 

Dept. 73-31, Bldg. 90, Plant A-l 
Burbank, CA 91520 


82 



DISTRIBUTION LIST (continued) 


Lockheed Palo Alto Research 
Laboratories 
Attn: D. Bushnell 
Palo Alto, CA 94304 

Lockheed Palo Alto Research 
Laboratories 
Attn: R. F. Hurtung 
Palo Alto, CA 94304 

Lockheed Missiles and Space Company 
Huntsville Research & Eng. Center 
Attn: W. Armstrong 
P. 0. Box 1103 
Huntsville, AL 35807 

MacNeal -Schwendler Corporation 
Attn: R. H. MacNeal 
7442 North Figueroa Street 
Los Angeles, CA 90041 

MARC Analysis Research Corporation 
Attn: P. V. Marcel 
260 Sheridan Avenue, Suite 314 
Palo Alto, CA 94306 

Northrop Space Laboratories 
Attn: Library 
3401 West Broadway 
Hawthorne, CA 90250 

North American Rockwell, Inc. 

Space & Information Systems Div. 
Attn: Library 
12214 Lakewood Blvd. 

Downey, CA 90241 

Norton Company 

Attn: Mr. George E. Buron 

Industrial Ceramics Div. 

Armore & Spectramic Products 
Worcester, MA 01606 

Norton Company 
Attn: Mr. Paul B. Gardner 
1 New Bond Street 
Industrial Ceramics Div. 

Worcester, MA 01606 


Republic Aviation 
Fairchild Hiller Corporation 
Attn: Library 

Farmington, Long Island, NY 

Rohr Industries 
Attn: Mr. John Meaney 
Foot of H Street 
Chula Vista, CA 92010 

Southwest Research Institute 

Attn: Dr. T. A. Cruse 

P. 0. Box 28510 

Culebra Road 

San Antonio, TX 78284 

Solar Turbine Inc. 

Attn: G. L. Padgett 
P. 0. Box 80966 
San Diego, CA 92138 

Cleveland State University 
Attn: Dr. P. Bel lino 
Department of Civil Engineering 
Cleveland, OH 44115 

Rockwell International Corporation 
Attn: J. Gausselin, D422/402AB71 
Los Angeles International Airport 
Los Angeles, CA 90009 

Rockwell International Corporation 

Rocketdyne Di vi si on 

Attn: J. F. Newell 

6633 Canoga Avenue 

Canoga Park, CA 91304 

Rockwell International Corporation 

Rocketdyne Division 

Attn: K. R. Rajagopal 

6633 Canoga Avenue 

Canoga Park, CA 91304 

TWA, Inc. 

Attn: Mr. John J. Morel li 
Kansas City International Airport 
P. 0. Box 20126 
Kansas City, M0 64195 


83 


DISTRIBUTION LIST (continued) 


Pratt & Whitney 

Attn: Library, MS 732-11 

P. 0. Box B269 

West Palm Beach, FL 33402 

Pratt & Whitney 

Attn: R. A. Marmol , MS 713-39 

P. 0. Box B269 

West Palm Beach, FL 33402 

Pratt & Whitney 

Attn: Library, MS 169-31 

400 Main Street 

East Hartford, CT 06108 

Pratt & Whitney 

Attn: E. S. Todd, MS 163-09 

400 Main Street 

East Hartford, CT 06108 

United Technologies Corporation 
Hamilton Standard 
Attn: Dr. R. A. Cornell 
Windsor Locks, CT 06096 

United Technologies Research Center 
Attn: Dr. A. Dennis, MS 18 
Silver Lane 

East Hartford, CT 06108 

Westinghouse R&D Center 
Attn: N. E. Dowling 
1310 Beulah Road 
Pittsburgh, PA 15235 

Westinghouse Adv. Energy Sys. Div. 

Attn: P. T. Falk 

Waltz Mill 

P. 0. Box 158 

Madison, PA 15663 

University of Akron 
Department of Civil Engineering 
Attn: Dr. Paul Chang 
Akron, OH 44325 


University of Arizona 
College of Engineering 
Attn: H. Kamel 
Tucson, AZ 87521 

University of California 
Department of Civil Engineering 
Attn: E. Wilson 
Berkeley, CA 94720 

University of Dayton 
Research Institute 
Attn: F. K. Bogner 
Dayton, OH 45409 

Georgia Institute of Technology 
School of Civil Engineering 
Attn: $. N. Atluri 
Atlanta, GA 30332 

Georgia Institute of Technology 
Attn: G. J. Simitsis 
225 North Avenue 
Atlanta, GA 30332 

Georgia Institute of Technology 
Attn: R. L. Cal son 
225 North Avenue 
Atl anta , GA 30332 

Lehigh University Institute of 
Fracture and Solid Mechanics 
Attn: G. T. McAllister 
Bethlehem, PA 18015 

Univ. of Illinois at Chicago Circle 
Dept, of Materials Engineering 
Attn: Dr. Robert L. Spilker 
Box 4348 

Chicago, IL 60680 

University of Kansas 
School of Engineering 
Attn: R. H. Dodds 
Lawrence , KS 66045 


84 



DISTRIBUTION LIST (continued) 


Massachusetts Institute of Technology 
Attn: J. Mar 
Cambridge, MA 02139 

State University of New York at Buffalo 
Department of Civil Engineering 
Attn: P. K. Banerjee 
Buffalo, NY 14225 

State University of New York at Buffalo 
Department of Civil Engineering 
Attn: R. Shaw 
Buffalo, NY 14214 

Purdue University 

School of Aeronautics & Astronautics 

Attn: C. T. Sun 

West Lafayette, IN 47907 

V. P. I. and State University 
Department of Engineering Mechanics 
Attn: R. H. Heller 
Blacksburg, VA 24061 


Rensselaer Polytechnic Institute 
Attn: R. Loewy 
Troy, NY 12181 

Stanford University 
School of Engineering 
Attn: T. J. R. Hughes 
Stanford, CA 94305 

Texas A&M University 
Aerospace Engineering Department 
Attn: W. E. Haisler 
College Station, TX 77843 

University of Texas at Austin 
College of Engineering 
Attn: J. T. Oden 
Austin, TX 7871 2-1 085 


85 


